© The Authors 2024

1 Introduction

The Yarkovsky-O’Keefe-Radzievskii-Paddack (YORP) effect is a thermal torque that can alter the spin state of an asteroid over time (Rubincam 2000; Vokrouhlickỳ & Čapek 2002; Vokrouhlicky et al. 2015), and is caused by the asymmetric re-emission of solar radiation by the irregular surface of the asteroid, resulting in a net torque that can spin up or spin down the asteroid’s rotation. The absorption of solar radiation makes no contributions to the YORP torque as it is averaged out over the spin and orbital periods for any asteroid shapes (Nesvornỳ & Vokrouhlickỳ 2008). So far, 11 asteroids showing time-varying rotational periods have been detected (Lowry et al. 2007; Taylor et al. 2007; Ďurech et al. 2022; Tian et al. 2022).

The YORP effect has important implications for an asteroid’s long-term rotational evolution. This effect can either spin down the asteroid to an extremely slow rotation, triggering a tumbling motion (Pravec et al. 2005), or spin up the asteroid to its rotation limit (e.g. rotational period of 2.2 h for rubble pile asteroids), leading to resurfacing (Sánchez & Scheeres 2020) and rotation disruption (Scheeres 2007; Fatka et al. 2020; Veras & Scheeres 2020). YORP-induced rotational disruption is supported by the observed asteroid pairs (Vokrouhlickỳ & Nesvornỳ 2008; Polishook 2014) and binary asteroids (Jacobson & Scheeres 2011; Delbo et al. 2011; Jacobson et al. 2013, 2016), including contact-binary asteroids (Rozek et al. 2019; Zegmott et al. 2021) and binary comets (Agarwal et al. 2020), which evolve under tidal effects and the binary YORP (BYORP) effect after the binary system is formed (Ćuk & Burns 2005; Steinberg & Sari 2011). Further potential observational evidence is the abnormal spin distribution and the obliquity distribution of near-Earth asteroids (Vokrouhlickỳ et al. 2003; Pravec et al. 2008; Rozitis & Green 2013b; Lupishko & Tielieusova 2014) and main belt asteroids (Lupishko et al. 2019), although a recent study points out that collisions might reproduce the observed distribution without the involvement of the YORP effect (Holsapple 2022). The YORP effect can influence the orbital evolution through the Yarkovsky effect, which is a thermal force that depends on the rotational state of the asteroid (Vokrouhlickỳ et al. 2000; Bottke et al. 2006). Therefore, understanding the YORP effect is important for correctly estimating the ages of asteroid families based on how much time is needed for the Yarkovsky effect to cause the observed orbital dispersion of their members from the original orbit following the disruption of their parent bodies (Vokrouhlickỳ et al. 2006; Ćuk et al. 2015; Carruba et al. 2014, 2015, 2016; Lowry et al. 2020; Marzari et al. 2020).

However, accurately calculating the YORP effect on a real asteroid remains a challenge, as it has been demonstrated to be highly sensitive to surface topology (Statler 2009; Breiter et al. 2009), such as uniform small-scale roughness (Rozitis & Green 2012), boulders (Golubov & Krugly 2012; Golubov et al. 2014, 2021; Ševeček et al. 2015; Golubov 2017; Golubov & Scheeres 2019; Golubov & Lipatova 2022), and craters (Zhou et al. 2022). Although the YORP torque caused by boulders and the tangential radiative force has been well studied, the YORP effect caused by concave structures has not yet been fully explored. More specifically, there are two kinds of concave structures on asteroid surfaces. The first one corresponds to craters, which result from impacts that an asteroid’s surface undergoes during its history and which can cover a large size range and be distributed in various ways. The second corresponds to surface roughness, which corresponds to uniformly distributed small-scale concave features that originate from the continuous effect of various processes that take place at the surface, such as thermal fatigue and space weathering. In our study, we consider both craters and surface roughness. While a pioneering study by Rozitis & Green (2012) used numerical simulations to investigate the effect of roughness, the high computational expense of such simulations prevents a comprehensive exploration of the functional dependence of this effect and its application to the global spin evolution of asteroid populations.

In addition to the precise calculation of the complete YORP effect, the long-term evolution of asteroids needs to account for stochastic collisions that affect this evolution, because each collision introduces a YORP torque due to the resulting crater. Bottke et al. (2015) performed a first study of the YORP effect accounting for craters caused by collisions, finding strong implications in the age estimate of asteroid families. However, their introduction of the concept of the stochastic YORP effect due to collisions assumed an arbitrary reset timescale for the YORP cycle. In reality, this reset timescale depends on the actual occurrence of each impact causing a crater on the asteroid’s surface and the resulting change in the YORP effect. In summary, collisions and the YORP and Yarkovsky effects are all coupled with each other in a way that is so far not well understood.

To account explicitly for the YORP torque caused by roughness and craters, Zhou et al. (2022) developed a semi-analytical model that is computationally efficient for the crater-induced YORP (CYORP) torque. The CYORP torque is defined as the torque difference between the crater and the flat ground: $$T_{\text{CYORP}}=T_{\text{crater}}-T_{\text{ground}}\text{.}$$(1)

In general, it takes the form of the following scaling rule with the radius of the crater R₀ and of the asteroid radius R_ast: $$T_{\text{CYORP}}=W\frac{\text{Φ}}{c}{R}_0^2{R}_{\text{ast}},$$(2)

where Φ is the flux of solar radiation at the asteroid’s semimajor axis, and W is the CYORP coefficient, which in turn depends on the obliquity and irregularity of the asteroid, and the depth-to-diameter ratio, location, thermal inertia, and albedo of the crater. As a first step, this model assumed a zero thermal conductivity.

Zhou et al. (2022) found that roughness or craters that cover 10% of the asteroid’s surface area could produce a CYORP torque comparable to the normal YORP (NYORP) torque, which arises from the macroscopic shape, with ignorance of the fine surface structure. Based on this number, which assumes that all craters have a depth-to-diameter ratio of 0.16, the reset of the YORP torque by the CYORP torque could be as short as 0.4 Myr.

However, the effects of finite thermal conductivity, self-radiation, and self-scattering were not considered in Zhou et al. (2022). In the present paper, we propose a semi-analytical model that accounts for the effects of self-sheltering, self-radiation, self-scattering, and non-zero thermal conductivity. This new model allows a more general exploration of the functional dependence of the CYORP effect. Moreover, as this semi-analytical model is much more efficient than a purely numerical one, it can be used to study the combined influence of the YORP effect, collisions, and the Yarkovsky effect by incorporating the CYORP effect into the standard evolution model of asteroid families. We assume that the craters considered in this work are on a convex asteroid. It is possible to apply the model derived in this work to a moderately concave asteroid by approximating concave structures as craters, but this is beyond the scope of this paper.

In this paper, we describe our analytical model for the temperature field in a crater in Sect. 2. In Sect. 3, we introduce the numerical model that we developed to validate the analytical model. The main results are shown in Sects. 4 and 5. In Sect. 4, we analyse the effects of self-sheltering, self-radiation, and self-scattering of the crater, and in Sect. 5 we show the results of the calculation of the CYORP torque as a function of different variables. For the purpose of illustration, in Sect. 6 we provide an example of the analysis of CYORP considering a specific real asteroid shape and its surface roughness as well as the consequential rotational evolution. In Sect. 7, we summarise our main findings and draw conclusions.

2 Analytical model

2.1 Calculation of temperature distribution in a crater

2.1.1 Linearized analytical solution

The temperature u for the surface and the layer beneath is governed by $$\frac{\partial u}{\partial t}=\frac{K}{C\rho }\frac{{\partial }^{2}u}{\partial z^2},$$(3)

with two boundary conditions, $${K\frac{\partial u}{\partial z}|}_{z=0}=F(t)-{e\sigma u^4|}_{z=0},$$(4) $${K\frac{\partial u}{\partial z}|}_{z\to \infty }=0,$$(5)

and a periodic condition, $${u|}_{t=2\pi /\omega }={u|}_{t=0},$$(6)

where t is the time, z is the depth of the crater, ω is the angular velocity, K is the thermal conductivity, C is the specific heat capacity, ρ is the bulk density of the asteroid, e is the emissivity, and σ is the Stefan-Boltzmann constant. In the following, we use the spin angle β = ωt to replace time for simplicity. The effect of the seasonal wave is marginal and is ignored in this work.

While the one-dimensional heat diffusion equation with the boundary condition of radiation has no complete analytical solution, it could be solved by linearizing the fourth power of the temperature, with the assumption that the temperature does not change significantly during a rotational period. For an arbitrary point in the hemispherical crater surface, the location of which is described by the polar angle θ and the azimuthal angle ϕ, the solution of the temperature is the real part of the expression $$u(\theta ,\varphi ,\beta )=u_0(\theta ,\varphi )+{\displaystyle \sum _{n=1}^{\infty }{u}_n}(\theta ,\varphi ){\text{e}}^{in\beta },$$(7)

with $$u_0(\theta ,\varphi )={\left(\frac{F_0(\theta ,\varphi )}{e\sigma }\right)}^{1/4},$$(8) $$u_n(\theta ,\varphi )=\frac{F_n(\theta ,\varphi )}{4e\sigma u_0{(\theta ,\varphi )}^3}\frac{{\text{e}}^{i{J}_n}}{\sqrt{2{\Theta }_n^2+2{\Theta }_n+1}},$$(9)

where $$l_n=\sqrt{\frac{n\omega }{2M}},$$(10) $${\text{Θ}}_n=\frac{K{l}_n}{4e\sigma u_0^3},$$(11) $$\tan J_n=-\frac{{\text{Θ}}_n}{{\text{Θ}}_n+1},$$(12)

with M = K/ρC. The function F_n is the nth-order of the Fourier series of the absorbed radiation flux F, which is expressed as $$F(\theta ,\varphi )={\displaystyle \sum _{n=0}^{\infty }{F}_n}(\theta ,\varphi ){\text{e}}^{in\beta }.$$(13)

We see that the only unknown variable is the absorbed radiation flux F. The absorbed radiation on a surface element contains three parts: solar radiation E(θ, ϕ, β), radiation from the crater itself H(θ, ϕ, β), and the scattering flux from the crater G(θ, ϕ, β), which leads us to $$\begin{array}{c}F(\theta ,\varphi )=(1-A)(E(\theta ,\varphi )+H(\theta ,\varphi )+G(\theta ,\varphi ))\\ =(1-A){\displaystyle \sum _{n=0}^{\infty }\left(E_n(\theta ,\varphi )+H_n(\theta ,\varphi )+G_n(\theta ,\varphi )\right)}{\text{e}}^{in\beta },\end{array}$$(14)

where E_n(θ, ϕ, β) and H_n(θ, ϕ, β) denote the nth-order Fourier modes of E(θ, ϕ, t) and H(θ, ϕ, t), respectively. Here the albedo A is assumed to be 0.1 for these three flux components for the sake of simplicity, although the albedo at the thermal-infrared wavelengths is almost zero. Following Eqs. (8), (9), and (14), we obtain $$u_0(\theta ,\varphi )={\left(\frac{(1-A)\left(E_0(\theta ,\varphi )+H_0(\theta ,\varphi )+G_0(\theta ,\varphi )\right)}{e\sigma }\right)}^{1/4},$$(15) $$u_n(\theta ,\varphi )=\frac{(1-A)\left(E_n(\theta ,\varphi )+H_n(\theta ,\varphi )+G_n(\theta ,\varphi )\right)}{4e\sigma u_0{(\theta ,\varphi )}^3}\frac{{\text{e}}^{i{J}_n}}{\sqrt{2{\text{Θ}}_n^2+2{\text{Θ}}_n+1}}.$$(16)

Therefore, to solve the temperature u(θ, ϕ, β), we need to obtain the Fourier series of E(θ, ϕ, β), H(θ, ϕ, β), and G(θ, ϕ, β).

Fig. 1 Three coordinate systems in this paper: coordinate system oxyz for calculating the illuminated domain in the crater, PABC for calculating the effective radiative force of an arbitrary surface element, and OXYZ for averaging the YORP torque over the spin and orbital motion.

2.1.2 Coordinate systems

We use three coordinate systems to calculate the radiation and the force received by the crater, as shown in Fig. 1. The coordinate system OXYZ is an inertial frame used to calculate the averaged YORP torque over the spin and orbital motion. Based on the axis OZ, we define the coordinate system oxyz fixed with the crater to calculate the instant solar flux. Finally, based on the axis oz, we define the coordinate system PABC to calculate the self-sheltering effect of the crater. The self-sheltering effect for a point in the crater refers to the non-working moment of the photons reabsorbed by the shelter (i.e. the crater itself), which leads to the effective radiation force of the surface being tilted relative to the surface normal with a modified magnitude. This self-sheltering effect on the point P depends on the geometry of the surrounding shelter. We adopt a simple sphere with a radius of R₁ to depict the crater shape.

In the coordinate system OXYZ, OZ points along the rotation axis and OXY is the equatorial plane. The axis OY is chosen such that the normal vector of the orbit plane lies in the plane OYZ. The axis OX follows the right-hand rule. There are three crucial vectors determining the CYORP torque expressed in polar coordinates in the coordinate system OXYZ: the crater position vector r₀ (denoted by α and β), the crater normal vector n₀ (described by α, β, δ, and Δ), and the solar position s (described bye and η), as shown in the right panel of Fig. 1. As the polar angle is between 0° and 180° by definition, we limit 0° < α < 180° and 0° < α−Δ < 180° in our results.

The coordinate system oxyz with the origin o located at the sphere centre is fixed with the crater in order to simplify the calculation of the solar flux on the crater. The axis oz points in the opposite direction to the surface normal vector n₀, which is also the symmetry axis of the spherical crater. The direction of axis oy is along e_OZ × e_oz and ox follows the right-hand rule. In this coordinate system, a crater with a depth of h can be defined as $$\text{𝒵}:=\left\{(x,y,z)\in {ℝ}^3\mid x^2+y^2+z^2=R_1,z\ge R_1\cos {\theta }_0\right\},$$(17)

with cos θ₀ = (R₁ − h)/R₁.

The coordinate system PABC with the origin P at a chosen point on the crater is used to calculate the self-sheltering effect. The axis PC is along the direction of Po, PA points in the direction of e_oz × e_PC, and PB follows the right-hand rule. In our code, the effective force felt by the point P is calculated first in the coordinate system PABC for simplicity and this is then transformed to the coordinate system oxyz by a rotation matrix.

2.1.3 Solar radiation E

In this section, we show how we derive the solar radiation received by an arbitrary point P in the crater, whose coordinates in oxyz system are r_P = (sin θ cos ϕ, sin θ sin ϕ, cos θ). The unit position vector of the Sun in the coordinate system OXYZ is described by $$s_{OXYZ}=(\cos \eta ,\cos ϵ\sin \eta ,\sin ϵ\sin \eta ),$$(18)

where ϵ is the obliquity and η is the angle of orbital motion. The transform matrix between the coordinate systems oxyz and OXYZ is set to $$ℛ=\left(\begin{array}{ccc}\cos \alpha \cos \beta & \cos \alpha \sin \beta & -\sin \alpha \\ \sin \beta & -\cos \beta & 0\\ -\cos \beta \sin \alpha & -\sin \alpha \sin \beta & -\cos \alpha \end{array}\right).$$(19)

Therefore, the coordinates of the vector s in the coordinate system oxyz is $$s_{\text{oxyz }}=\left(\begin{array}{c}-\sin \alpha \sin ϵ\sin \eta +\cos \alpha \cos \beta \cos \eta +\cos \alpha \sin \beta \cos ϵ\sin \eta \\ \sin \beta \cos \eta -\cos \beta \cos ϵ\sin \eta \\ -\sin \alpha \cos \beta \cos \eta -\sin \alpha \sin \beta \cos ϵ\sin \eta -\cos \alpha \sin ϵ\sin \eta \end{array}\right).$$(20)

On the other hand, we can use the angle λ and ϕ′ to represent s_oxyz: $$s_{\text{oxyz }}=\left(\sin \lambda \cos {\varphi }^{\prime },\sin \lambda \sin {\varphi }^{\prime },-\cos \lambda \right),$$(21)

such that $$\cos \lambda =\sin \alpha \cos \beta \cos \eta +\sin \alpha \sin \beta \cos ϵ\sin \eta +\cos \alpha \sin ϵ\sin \eta ,$$(22)

and $$\tan {\varphi }^{\prime }=\frac{\sin \alpha \sin ϵ\sin \eta -\cos \alpha \cos \beta \cos \eta -\cos \alpha \sin \beta \cos ϵ\sin \eta }{\sin \beta \cos \eta +\cos \beta \cos ϵ\sin \eta }.$$(23)

The absorbed radiation flux is $$\begin{array}{ll}E(\theta ,\varphi )\hfill & =(1-A)\text{Φ}H(\cos \lambda )H(w)\left(-r_P\cdot s\right)\hfill \\ \hfill & =(1-A)\text{Φ}H(\cos \lambda )H(w)\cdot \left[m_1\cos \left(\beta -{\beta }_1\right)+m_2\right]\hfill \end{array}$$(24)

where $$w=\cos 2\lambda \cos \theta +\sin {\theta }_0-\sin 2\lambda \sin \theta \cos \left(\varphi -{\varphi }^{\prime }\right),$$(25) $$\tan {\beta }_1=\frac{(\sin \alpha \cos \theta -\cos \alpha \sin \theta \cos \varphi )\cos ϵ\sin \eta -\sin \theta \sin \varphi \cos \eta }{(\sin \alpha \cos \theta -\cos \alpha \sin \theta \cos \varphi )\cos \eta +\cos ϵ\sin \eta \sin \theta \sin \varphi },$$(26) $$\begin{array}{l}{m}_1={\left({\cos }^2\eta +{\cos }^2ϵ{\sin }^2\eta \right)}^{1/2}\\ \cdot ({\cos }^2\theta {\sin }^2\alpha -\frac{\cos \varphi \sin 2\alpha \sin 2\theta }{2}\\ {+\left({\cos }^2\alpha {\cos }^2\varphi +{\sin }^2\varphi \right){\sin }^2\theta )}^{1/2},\end{array}$$(27) $$m_2=\cos \alpha \cos \theta \sin ϵ\sin \eta +\cos \varphi \sin \alpha \sin ϵ\sin \eta \sin \theta \text{. }$$(28)

Here, H is the Heaviside function defined by $$H(x):=\{\begin{array}{l}1,x>0\hfill \\ 0,\text{ else.}\hfill \end{array}$$(29)

Using the substitution β′ = β − β₁, the absorbed radiation flux has the form $$E(\theta ,\varphi )=\{\begin{array}{l}(1-A)\text{Φ}\left(m_1\cos {\beta }^{\prime }+m_2\right),{{\beta }^{\prime }}_{\min }<\beta <{{\beta }^{\prime }}_{\max }\hfill \\ 0,\text{ else. }\hfill \end{array}$$(30)

Expanding Eq. (30) in Fourier series, we obtain $$E(\theta ,\varphi )=E_0+\frac{(1-A)\text{Φ}}{\pi }{\displaystyle \sum _{n=1}^{\infty }\left[C_n\cos \left(n{\beta }^{\prime }\right)+S_n\sin \left(n{\beta }^{\prime }\right)\right]},$$(31)

with $$E_0={\frac{(1-A)\text{Φ}}{2\pi }\left(m_1\sin \beta +m_2\beta \right)|}_{{{\beta }^{\prime }}_{\min }}^{{{\beta }^{\prime }}_{\max }},$$(32) $$\begin{array}{l}{C}_n={\displaystyle {\int }_{{\beta }^{\prime }}^{{{\beta }^{\prime }}_{\min }}\left(m_1\cos {{\beta }^{\prime }}_{\max }+m_2\right)}\cos \left(n{\beta }^{\prime }\right)\text{d}{\beta }^{\prime }\\ =\{\begin{array}{l}\left(\frac{m_1\sin 2{\beta }^{\prime }+2m_1{\beta }^{\prime }}{4}+m_2\sin \beta \right)\begin{array}{l}{{\beta }^{\prime }}_{\max }\hfill \\ {\beta }_{\min }\hfill \end{array},n=1\hfill \\ \left(\frac{m_{1}n\cos {\beta }^{\prime }\sin n{\beta }^{\prime }-m_1\cos n{\beta }^{\prime }\sin {\beta }^{\prime }}{n^2-1}+\frac{m_2\sin \left(n{\beta }^{\prime }\right)}{n}\right){{\beta }^{\prime }}_{\min },n>1,\hfill \end{array}\end{array}$$(33)

and $$\begin{array}{c}{S}_n={\displaystyle {\int }_{{\beta }^{\prime }}^{{\beta }_{{\min }^{\prime }}}\left(m_1\cos {\beta }_{{\max }^{\prime }}+m_2\right)}\sin \left(n{\beta }^{\prime }\right)\text{d}{\beta }^{\prime }\\ =\{\begin{array}{l}-\left(\frac{m_1\cos {\beta }^{\prime 2}}{2}+m_2\cos {\beta }^{\prime }\right)|\begin{array}{c}{\beta }_{\max }^{\prime }\\ {\beta }_{\min }^{\prime }\end{array},n=1\hfill \\ \left(\frac{m_{1}n\cos {\beta }^{\prime }\cos n{\beta }^{\prime }+m_1\sin {\beta }^{\prime }\sin n{\beta }^{\prime }}{1-n^2}-\frac{m_2\cos n{\beta }^{\prime }}{n}\right)|{\beta }_{\max }^{\prime },n>1.\hfill \end{array}\end{array}$$(34)

The nth coefficient of the Fourier series of E(θ, ϕ) is $$E_n=\sqrt{C_n^2+S_n^2}\cdot {\text{e}}^{i{\text{Φ}}_n},$$(35)

with $$\tan {\text{Φ}}_n=-\frac{S_n}{C_n}.$$(36)

2.1.4 Self-heating effect

Due to the concave structure of the crater, the surface element in the crateralso receives photons emitted by the crateritself, which is a process referred to as ‘self-heating’. In this section, we derive the self-radiation H and the self-scattering G as a function of the position in the crater.

Let us consider a surface element dS₁ receiving the radiation from another surface element dS₂. In the reference frame oxyz, the positions of dS₁ and dS₂ are $$r_1=R_1\left(\sin {\theta }_1\cos {\varphi }_1,\sin {\theta }_1\sin {\varphi }_1,\cos {\theta }_1\right),$$(37) $$r_2=R_1\left(\sin {\theta }_2\cos {\varphi }_2,\sin {\theta }_2\sin {\varphi }_2,\cos {\theta }_2\right)$$(38)

The displacement from dS₁ to dS₂ is $$r_{1,2}=r_2-r_1.$$(39)

The incident angle ζ₁ and the emission angle ζ₂ are defined as $$\cos {\zeta }_1=-r_1\cdot r_{1,2}/R_1{r}_{1,2}$$(40) $$\cos {\zeta }_2=r_2\cdot r_{1,2}/R_1{r}_{1,2},$$(41)

respectively. The radiation flux at the location of dS₁ produced by the thermal radiation of dS₂ is $$H_{1,2}=\frac{e\sigma u_2^4}{\pi }\frac{\cos {\zeta }_1\cos {\zeta }_2}{r_{1,2}^2}\text{d}{S}_2$$(42)

Substituting Eqs. (37)–(41) into Eq. (42), we obtain $$H_{1,2}=\frac{e\sigma u_2^4}{4\pi R_1^2}\text{d}{S}_2.$$(43)

For an arbitrary point, the radiation flux caused by the whole crater is $$H(\theta ,\varphi )=\underset{\mathcal{Z}}{{\displaystyle \int }}\frac{e\sigma u{\left({\theta }^{\prime },{\varphi }^{\prime }\right)}^4}{4\pi }\sin {\theta }^{\prime }\text{d}{\theta }^{\prime }\text{d}{\varphi }^{\prime },$$(44)

where 𝒵 is the crater surface, and is defined as $$\mathcal{Z}:=\left\{(x,y,z)\in {ℝ}^3\mid r=R_1,\theta \in \left(0,\pi /2-{\theta }_0\right),\varphi \in (0,2\pi )\right\}.$$(45)

Similarly, we obtain the self-scattering term: $$G(\theta ,\varphi )=\underset{\mathcal{Z}}{{\displaystyle \int }}\frac{AE\left({\theta }^{\prime },{\varphi }^{\prime }\right)}{4\pi }\sin {\theta }^{\prime }\text{d}{\theta }^{\prime }\text{d}{\varphi }^{\prime }.$$(46)

Therefore, both H(θ, ϕ) and G(θ, ϕ) can be expressed in terms of u(θ, ϕ) and E(θ, ϕ). As E(θ, ϕ) is derived in Sect. 2.1.3, the only unknown variable is the temperature u(θ, ϕ), the solution for which is discussed in the following section.

2.2 Solution for temperature

We obtained the Fourier series of the solar radiation flux (Sect. 2.1.3) and the radiation flux produced by the crater (Sect. 2.1.4), which allows us to return to Eq. (14) to solve the temperature distribution u(θ, ϕ). We note that the self-radiation term H contains the unknown temperature distribution that is to be solved.

The Fourier coefficients of the temperature of the crater follow $$\begin{array}{c}\sigma u_0^4(\theta ,\varphi )=(1-A)\left(E_0(\theta ,\varphi )+G_0(\theta ,\varphi )+H_0(\theta ,\varphi )\right)\\ \sigma u_n(\theta ,\varphi )u_0^3(\theta ,\varphi )\frac{2\sqrt{2{\text{Θ}}_n^2+2{\text{Θ}}_n+1}}{(1-A){\text{e}}^{i{J}_n}}=E_n(\theta ,\varphi )+G_n(\theta ,\varphi )+H_n(\theta ,\varphi ),\end{array}$$(47)

with $$\begin{array}{c}{G}_n(\theta ,\varphi )=\underset{\mathcal{Z}}{{\displaystyle \int }}\frac{A{E}_n\left({\theta }^{\prime },{\varphi }^{\prime }\right)}{4\pi }\sin {\theta }^{\prime }\text{d}{\theta }^{\prime }\text{d}{\varphi }^{\prime },\\ H_0(\theta ,\varphi )=\underset{\mathcal{Z}}{{\displaystyle \int }}\frac{e\sigma u_0{\left({\theta }^{\prime },{\varphi }^{\prime }\right)}^4}{4\pi }\sin {\theta }^{\prime }\text{d}{\theta }^{\prime }\text{d}{\varphi }^{\prime },\\ H_n(\theta ,\varphi )=\underset{\mathcal{Z}}{{\displaystyle \int }}\frac{e\sigma u_0{\left({\theta }^{\prime },{\varphi }^{\prime }\right)}^3{u}_n\left({\theta }^{\prime },{\varphi }^{\prime }\right)}{\pi }\sin {\theta }^{\prime }\text{d}{\theta }^{\prime }\text{d}{\varphi }^{\prime }.\end{array}$$(48)

Here, terms G and H are the scattering flux and self-radiation flux, respectively.

2.2.1 Solution of a general form

Equation (47) have a general form: $$f(\theta ,\varphi )=g(\theta ,\varphi )+C\underset{\mathcal{Z}}{{\displaystyle \int }}f\left({\theta }^{\prime },{\varphi }^{\prime }\right)\sin {\theta }^{\prime }\text{d}{\theta }^{\prime }\text{d}{\varphi }^{\prime }.$$(49)

By setting $$k=\underset{\mathcal{Z}}{{\displaystyle \int }}f\left({\theta }^{\prime },{\varphi }^{\prime }\right)\sin {\theta }^{\prime }\text{d}{\theta }^{\prime }\text{d}{\varphi }^{\prime },$$(50)

we have $$f(\theta ,\varphi )=g(\theta ,\varphi )+C\cdot k.$$(51)

Substituting Eq. (51) into Eq. (50), we obtain $$k=\underset{\mathcal{Z}}{{\displaystyle \int }}\left(g\left({\theta }^{\prime },{\varphi }^{\prime }\right)+C\cdot k\right)\sin {\theta }^{\prime }\text{d}{\theta }^{\prime }\text{d}{\varphi }^{\prime }.$$(52)

Rearranging Eq. (52), we find k is $$k=\frac{{{\displaystyle \int }}_{\mathcal{Z}}g\left({\theta }^{\prime },{\varphi }^{\prime }\right)\sin {\theta }^{\prime }\text{d}{\theta }^{\prime }\text{d}{\varphi }^{\prime }}{1-C{{\displaystyle \int }}_{\mathcal{Z}}\sin {\theta }^{\prime }\text{d}{\theta }^{\prime }\text{d}{\varphi }^{\prime }},$$(53)

with which f(θ, ϕ) is solved out by Eq. (51).

2.2.2 Solution for temperature

In the case of u₀, $$f(\theta ,\varphi )=e\sigma u_0^4/(1-A),$$(54) $$g(\theta ,\varphi )=E_0(\theta ,\varphi )+G_0(\theta ,\varphi ),$$(55) $$C=\frac{1}{4\pi },$$(56)

and in the case of u_n, $$f(\theta ,\varphi )=e\sigma u_0^3{u}_n\frac{2\sqrt{2{\text{Θ}}_n^2+2{\text{Θ}}_n+1}}{(1-A){\text{e}}^{i{J}_n}},$$(57) $$g(\theta ,\varphi )=E_n(\theta ,\varphi )+G_n(\theta ,\varphi ),$$(58) $$C=\frac{{\text{e}}^{i{J}_n}}{2\pi \sqrt{2{\text{Θ}}_n^2+2{\text{Θ}}_n+1}}.$$(59)

2.3 The CYORP torque

The CYORP torque is defined as the YORP torque difference between the crater and the flat ground: $$T_{\text{CYORP }}=T_{\text{crater }}-T_{\text{ground }}.$$(60)

2.3.1 The YORP torque of the crater

The average radiative torque produced by the crater should be calculated in the inertial frame OXYZ: $$T_{\text{crater }}=r_{0,OXYZ}\times f_{OXYZ}.$$(61)

Here, f_OXYZ is the radiative force, and $$f_{OXYZ}=\underset{\mathcal{W}}{{\displaystyle \int }}\frac{e\sigma T{(\theta ,\varphi )}^4}{c}{n}_f(\theta ,\varphi )\text{d}S,$$(62)

where n_f is the corrected force direction vector for each surface element. If there is no shelter around the surface element, n_f is equal to the surface unit normal vector n₀. However, the surface element in a crater has a sky sheltered by other elements, resulting in the reabsorption of the emitted photons along the direction of the shelter.

For the surface element dS(θ, ϕ), the radiative force is $$f=-\underset{\mathscr{H}}{{\displaystyle \int }}\frac{\text{Φ}}{\pi c}\cos {\theta }^{\prime }\left(\begin{array}{l}\sin {\theta }^{\prime }\cos {\varphi }^{\prime }\hfill \\ \sin {\theta }^{\prime }\sin {\varphi }^{\prime }\hfill \\ \cos {\theta }^{\prime }\hfill \end{array}\right)\sin {\theta }^{\prime }\text{d}{\theta }^{\prime }\text{d}{\varphi }^{\prime }.$$(63)

Without sheltering (e.g. for convex asteroids), ℋ is replaced by the hemisphere (i.e. θ ∈ (0,π/2), ϕ ∈ (0, 2π)). In this case, the force is reduced to 2Φn₀/3c.

2.3.2 The YORP torque of the flat portion of the surface

The absorbed radiation flux for a flat ground with the normal vector n₀ is $$\begin{array}{c}{E}_{\text{ground }}=(1-A)\text{Φ}H(\cos \lambda )\cos \lambda \\ =\{\begin{array}{c}(1-A)\text{Φ}\left(m_1\cos {\beta }^{\prime }+m_2\right),{\beta }_{\min }^{\prime }<{\beta }^{\prime }<{\beta }_{\max }^{\prime }\\ 0,\text{ else. }\end{array}\end{array}$$(64)

Here, β′ = β −β₁, with $$\tan {\beta }_1=\frac{\sin \alpha \cos ϵ\sin \eta }{\sin \alpha \cos \eta },$$(65) $$m_1={\left({\sin }^2\alpha \left({\cos }^2\eta +{\cos }^2ϵ{\sin }^2\eta \right)\right)}^{1/2},$$(66) $$m_2=\cos \alpha \sin ϵ\sin \eta ,$$(67)

and β_min and β_max are the negative and positive values of arccos(−m₂/m₁), respectively.

3 Numerical model for examination

In the above analytical method, we adopted the assumption of a “small” temperature variation during a rotation period, which allows us to linearize the biquadrate of the temperature (i.e. $u^4~u_0^3{{\displaystyle \Sigma }}^{\text{}}{u}_n{\text{e}}^{\beta n\text{i}}$). This assumption is equivalent to a high value of the thermal parameter, which is defined as $$\text{Γ}=\frac{\sqrt{C\rho \omega \kappa }}{{(ϵ\sigma )}^{1/4}{(1-A)}^{3/4}{\text{Φ}}^{3/4}}.$$(68)

In the case of a low thermal parameter, the analytical model should be used with caution. To examine the appropriate range of the thermal inertia for which our analytical model is valid, we built a one-dimensional thermophysical numerical model to perform cross-validation.

3.1 Numerical model

We used a finite difference numerical method to solve Eq. (3) with the second-order Crank-Nicholson scheme: $$\begin{array}{r}\hfill c_{\text{n}}{u}_{i,j+1,k+1}-\left(2c_{\text{n}}+1\right)u_{i,j,k+1}+c_{\text{n}}{u}_{i,j-1,k+1}=\\ \hfill -c_{\text{n}}{u}_{i,j+1,k}+\left(2c_{\text{n}}-1\right)u_{i,j,k}-c_{\text{n}}{u}_{i,j-1,k}.\end{array}$$(69)

Here, u_i,j,k represents the temperature at the depth of (j − 1)δz below the ith facet at the kth time step, where i, j, and k are integrals starting from 1 to i_max, j_max, and k_max. The coefficient c_n is $$c_{\text{n}}=\frac{a\delta t}{2{(\delta x)}^2}.$$(70)

The value of c_n should be smaller than 0.5 for the convergence of the iteration.

The surface temperature is determined by the first boundary condition (Eq. (4)): $$(1-A)\left(E_{i,k}+H_{i,k}+S_{i,k}\right)-\sigma u_{i,1,k+1}^4=K\frac{u_{i,1,k+1}-u_{i,2,k}}{\delta z},$$(71)

which can be solved by a Newtonian-Raphson iterative method. The solar flux E_i,k on the ith facet at the kth time step is $$E_{i,k}=\{\begin{array}{l}{\text{Φ}}_0{n}_i\cdot s_{k}H\left(n_i\cdot s_k\right),\text{ unsheltered }\hfill \\ 0,\text{ sheltered, }\hfill \end{array}$$(72)

where Φ₀ is 1364 W m⁻² at the distance of 1 au to the Sun. Here, n_i is the normal vector of the ith facet and s_k is the unit solar position vector at the kth time step. Whether or not the facet is sheltered is judged according to the projections of other facets on the plane of the ith facet along the solar position vector s_k. The self-radiation term H_i,k is the sum of radiation from other facets: $$H_{i,k}=\underset{i^{\prime }\ne i}{{{\displaystyle \Sigma }}^{\text{}}}e\sigma u_{i^{\prime },1,k}^4\frac{-\left(n_i\cdot r_{i,i^{\prime }}\right)\left(n_{i^{\prime }}\cdot r_{i,i^{\prime }}\right)}{\pi r_{i,i^{\prime }}^2}{S}_{i^{\prime }}$$(73)

where r_i,i′, is the vector from the centre of the ith facet to the centre of i′th facet, and S_i′ is the area of the ith facet. The scattering term S_i,k is given by $$S_{i,k}=\underset{i^{\prime }\ne i}{{{\displaystyle \Sigma }}^{\text{}}}A{E}_{i^{\prime }}\frac{-\left(n_i\cdot r_{i,i^{\prime }}\right)\left(n_{i^{\prime }}\cdot r_{i,i^{\prime }}\right)}{\pi r_{i,i^{\prime }}^2}{S}_{i^{\prime }}.$$(74)

The second boundary condition (Eq. (5)) translates into $$u_{i,j_{\max },k+1}=u_{i,j_{\max }-1,k+1}.$$(75)

Combining Eqs. (69), (71), and (75), we can obtain the solution for the temperature in the next time step u_{i,j,k + 1}·,based on the temperature in the current time step u_i,j,k We adopted an initial temperature of u_i,j,0 = 280 K. The maximum depth is set to be a few thermal skin depths $\sqrt{K/C\rho \omega }$ and the total number of layers is set as j_max = 50. In order to make sure that the surface temperature is in an equilibrium state, we set the total time as k_maxδt ~ 20 spin periods. We divided the considered crater into about 1000 facets and solve the temperature for each facet using the above method.

Fig. 2 Radiative force of the total crater averaged over a rotational period (8 h by default) obtained from the analytical method (solid lines) and the numerical method (dashed lines), as a function of the crater colatitude α for K = 2.65 Wm−1 K−1 (left panel) and K = 0.1 Wm−1 K−1 (right panel). The x, y, and 𝓏 components of the radiative force are shown in blue, red, and green, respectively. The distance of the crater here is 1 au from the Sun.

3.2 Comparison with the analytical model

The thermal parameters of asteroids can vary widely depending on their composition and structure. For example, the thermal conductivity of a porous material is much lower than that of a dense metal. The thermal conductivity of stony asteroids, which are made mostly of silicates, can range from about 0.1 Wm K⁻¹ to 1 WmK⁻¹, while the thermal conductivity of metallic asteroids is generally much higher, in the range of 20–50 WmK⁻¹. Asteroids that are composed of a mixture of rock and metal will have thermal conductivity values between those of pure rock and pure metal.

Here, we test three typical types of asteroid materials: regolith, solid basalt, and metal, whose properties are shown in Table 1. We calculate the radiative force of the total crater averaged over a rotational period (8 h by default), as a function of the colatitude of the crater. The craters in the test are placed at 1 au from the Sun. For simplicity, we set the obliquity to ϵ = 0. The results computed from the analytical method and the numerical method are shown in Fig. 2. We can see that the analytical result is consistent with the numerical result to a high degree, while a large discrepancy shows up when the thermal conductivity decreases to 0.1 Wm⁻¹ K⁻¹. This coincides with our prerequisite of the application of our analytical method, which is that the temperature variation should be small. Therefore, our method is appropriate for solid basalt and metal materials.

Regarding regolith material, with a thermal conductivity of as low as 0.001 Wm⁻¹ K⁻¹, we test the model described in Zhou et al. (2022), where zero thermal conductivity is assumed. The result is shown in Fig. 3, which indicates that this latter model works well for regolith materials. Therefore, for three materials representing asteroid surfaces, our two methods, namely the one in the present work and that described in Zhou et al. (2022), behave well in modelling the YORP effect.

4 Discussion on self-modification effects

For a concave structure, there are three self-modification effects: the self-sheltering, self-radiation, and self-scattering effect. The first one refers to the radiative force modification on the surface element due to the re-absorption of photons by the crater. The second and third ones refer to the temperature increase due to the emitted photons and scattered photons from the crater itself. Previous research shows that these self-modification effects could be important for the YORP torque of the crater (Statler 2009; Rozitis & Green 2012, 2013a), but a quantitative description is still lacking. For example, it is not clear how deep the crater needs to be so that these self-modification effects can no longer be ignored. This is crucial to the validity of the commonly made assumption that the asteroid can be considered as a moderately convex shape with the craters being overlooked when evaluating the YORP torque. An essential metric in this respect is the depth-to-diameter ratio of the crater. We explore the radiative force (leading to the YORP torque) of a crater based on its depth-to-diameter ratio under varying conditions - one for each self-modification effect. Figure 4 shows our findings, situating the craters at the equator (α = π/2) with a spin axis of zero-degree obliquity.

Our analysis reveals that both self-radiation and self-scattering amplify the force, raising the temperature by as much as 50% when the crater mirrors a hemisphere. In contrast, self-sheltering diminishes the force by up to 30%. Notably, the impact of self-scattering remains negligible for typical asteroid surface albedos ranging between 0.1 and 0.2. When integrating all self-modification effects, the radiative force increases by 15%.

Our result also shows that when the depth-to-diameter ratio h/D₀ < 0.05, the force increase is less than 1%. Therefore, for those shallow concave structures with h/D₀ < 0.05, no self-modification effects are needed in the YORP model for it to remain accurate. These can then be efficiently approximated as flat surfaces.

Fig. 3 Same as Fig. 2 but for K = 0 Wm−1 K−1 in the analytical model and K = 0.001 Wm−1 K−1 in the numerical model.

Fig. 4 Ratio of the radiative force including self-modification effects to that without self-modification effects, as a function of crater depth-to-diameter ratio and accounting for different self-modification effects.

5 Analysis of the CYORP torque

As shown in a previous work (Zhou et al. 2022), the CYORP torque depends on many factors, including the depth-to-diameter ratio and location (or the colatitude α for example) of the crater as well as the obliquity and thermal parameter of the asteroid. In this section, we discuss the dependence of the CYORP torque on these factors. In the following, except for in Sect. 5.1, we assume the depth-to-diameter ratio to be ~0.16.

5.1 Depth-to-diameter ratio

Asteroid craters exhibit a range of distinct features in size and shape, with diameters ranging from a few centimetres to several kilometres for large asteroids. These craters generally display bowl-shaped structures, containing central peaks and terraced walls when produced in the gravity regime. To simplify the modelling, a semi-sphere approximation is often used to represent the shape of craters.

According to the definition of CYORP torque, if the depth-to-diameter ratio reaches zero, the CYORP torque is zero (Eq. (60)). Figure 5 shows the CYORP torques generated by craters with various depth-to-diameter ratios. The parameters δ = ∆ = π/6, K = 2.65 Wm⁻¹ K⁻¹, and α = π/2 are used. Higher depth-to-diameter ratios correspond to larger spin components of the CYORP torque. A crater with h/D₀ < 0.05 produces an insignificant CYORP torque, which may be disregarded. Furthermore, the depth-to-diameter ratio also impacts the obliquity component, influencing both the torque magnitude and shape of the torque curve. For instance, when the depth-to-diameter ratio is low, the asymptotic obliquity is 90°, while for higher depth-to-diameter ratios, new asymptotic obliquities arise around 0° and 180°.

5.2 Crater latitude α and asteroid obliquity ϵ

Figure 6 displays the CYORP torque components as a function of the crater latitude and asteroid obliquity. The values of the parameters δ and ∆ are set to a representative value of π/4. The spin component W_𝓏 exhibits symmetry about the axis 𝓏 = 90°, while the obliquity component W_𝓏 is anti-symmetric. The minimum and maximum values of the torque occur when the obliquity is 0^° or 90^°, with the absolute value of these extrema reaching up to 0.02. The coefficient of the obliquity component of the CYORP torque is considerably smaller, with a maximum value of 0.004.

For comparison, the typical value of the normal YORP spin coefficient is 0.005 for type I/II and <0.001 for type III/IV aster-oids¹. The ratio of the CYORP torque to the normal YORP torque scales as $$\frac{T_{\text{CYORP}}}{T_{\text{NYORP}}}=\frac{W_{\text{CYORP}}{S}_{\text{crater }}}{W_{\text{YORP}}{S}_{\text{asteroid }}}.$$(76)

Setting the ratio to 1, we find that the total area of concave structures needs to be as large as 1/4 and 1/20 of the asteroid surface area for type I/II and type III/IV asteroids, respectively.

Fig. 5 Spin component (left panel) and the obliquity component (right panel) of the CYORP coefficient as a function of the asteroid obliquity e, accounting for different depth-to-diameter ratios, which are denoted by different colours.

Fig. 6 Spin component (left panel) and the obliquity component (right panel) of the CYORP coefficient as a function of the asteroid obliquity e and the crater colatitude α(> ∆), for K = 2.65Wm−1 K−1. Here, d/D0 = 0.16 and δ = ∆ = π/4.

5.3 Thermal parameter

When the asteroid rotates quickly and has high heat conductivity, a higher value of the thermal parameter arises, resulting in a smaller variation in temperature. To explore the role of the thermal parameter, we employ the same parameters as in Sect. 5.2, but with K = 40 Wm⁻¹ K⁻¹ for metal materials, and the resulting CYORP torques are depicted in Fig. 7. The comparison with Fig. 6 reveals that the spin component remains relatively unchanged, while the obliquity component displays significant variation. This observation aligns with the prior assertion that the thermal parameter mainly impacts the obliquity component (Vokrouhlickỳ & Čapek 2002). In the regime of high thermal conductivity, the obliquity component diminishes as the thermal conductivity increases.

5.4 Irregularity δ and ∆

The angular parameters δ and ∆ are used to describe the irregularity of the asteroid, where δ = 0 and ∆ = 0 correspond to a perfect sphere. We explore the CYORP torque as a function of δ and ∆ with fixed asteroid obliquity and crater colatitude of π/4. The results are presented in Fig. 8. We can see that for the spin component, δ controls the torque magnitude while ∆ controls the torque direction.

Zhou et al. (2022) demonstrates that the CYORP torque vanishes for δ = 0. However, in the presence of finite thermal inertia, the obliquity component of the CYORP torque arises while the spin component remains negligible. Figure 9 illustrates the variation of the CYORP obliquity component with the asteroid obliquity and the crater colatitude when δ and ∆ are both zero. The CYORP torque has a tendency to lead the asteroid obliquity to 90° when the crater is near the poles, while it leads to an asymptotic obliquity of 0° or 180° when the crater is near the equator.

Fig. 7 Same as Fig. 6 but for K = 40 Wm−1 K−1.

Fig. 8 Spin component (left panel) and the obliquity component (right panel) of the CYORP coefficient as a function of ∆(<α) and δ, for K = 2.65 Wm−1 K−1. Here, d/D0 = 0.16 and a = ϵ = π/6.

6 Application of the CYORP effect on a real asteroid

6.1 Roughness

The surface roughness of asteroids is produced by several processes, including micrometeorite impacts, thermal fatigue, ejecta, or outgassing. It was found that the YORP torque is extremely sensitive to the small-scale surface structures (Statler 2009; Breiter et al. 2009). The microscopic beaming effect of regolith grain-size-scale roughness (<1 mm) was shown – using the Hapke reflectance and emissivity model (Breiter & Vokrouhlickỳ 2011) – to have a marginal influence on the YORP effect. The transverse heat conduction across thermal skin depth (~1 cm) causes an asymmetric thermal emission of the east and west sides of a boulder, giving rise to a systematic positive YORP torque (Golubov & Krugly 2012; Golubov & Lipatova 2022). The importance of the various self-modification effects of a concave feature on the surface was considered gradually and numerical approaches were taken to study it (Statler 2009; Rozitis & Green 2012, 2013a). It was found that the concave feature of surface roughness could dampen the YORP torque by tens of percent. While the pioneering studies by Rozitis & Green (2012) shed light on the effects of roughness, the computational expense and difficulty in studying the functional dependence means that there are severe limitations to the numerical method.

In contrast, the analytical method that we introduce in the present study, and its computational efficiency, simplify the application of roughness-induced YORP effects to real asteroids or binary systems. Given the objective of our semi-analytical method to provide a rapid assessment of the impact of surface roughness, it is particularly well-suited for models of asteroids with rough convex shapes derived from light-curve observations. However, when dealing with high-resolution shape models, especially those that possess a resolution of a few centimetres (the scale of the thermal skin), a 3D thermal model becomes essential for accurate calculations, owing to the presence of the tangential YORP (TYORP) effect, which accounts for the transverse heat conduction inside a boulder.

For illustrative purposes, to demonstrate the application of our method, we randomly selected the main belt asteroid (272) Antonia as an example. This choice is representative of the majority of asteroids lacking detailed information obtained through in situ observations. We used the shape model obtained from Hanus et al. (2013). To optimise the performance of our model, we assume a relatively high thermal conductivity of 1 Wm⁻¹ K⁻¹. We uniformly distributed the roughness across the asteroid’s surface and investigated the total YORP torque (NYORP + CYORP). To do so, we input the normal vector and position vector of surface elements in the shape model to our CYORP model to obtain the CYORP torque coefficient of each surface element. We then used the area of each surface element to calculate the CYORP torque and sum up all CYORP torques generated by all surface elements. The depth-to-diameter is assumed to be 0.5, following the assumption made by Rozitis & Green (2012). This has been compared with the sole NYORP torque. The result is depicted in Fig. 10. Our findings confirm that the roughness-induced CYORP torque damps the normal YORP torque. Specifically, for the asteroid (272) Antonia, the spin component of the torque experiences a damping effect of approximately 35%, while the obliquity component is damped by approximately 64%.

Fig. 9 Obliquity component (right panel) of the CYORP coefficient as a function of the asteroid obliquity e and the crater colatitude α for K = 2.65 Wm−1 K−1. Here, d/D0 = 0.16 and δ = ∆ = 0.

6.2 Rotational evolution

The rotational dynamics of asteroids are primarily governed by two key processes: collisions and the YORP effect. The timescale for reorientation of the spin axis resulting from angular momentum transfer during a collision can be expressed as (Farinella et al. 1998) $${\tau }_{\text{imp,re-ori }}=746{\left(\frac{R_{\text{ast }}}{1 \text{km}}\right)}^{4/3}{\left(\frac{\omega }{3\times {10}^{-4} {\text{s}}^{-1}}\right)}^{5/6}\text{Myr}.$$(77)

On the other hand, the typical timescale associated with the YORP effect is approximately given by: $${\tau }_{\text{YORP}}~4{\left(\frac{R_{\text{ast}}}{1 \text{km}}\right)}^2\left(\frac{\omega }{3\times {10}^{-4} {\text{s}}^{-1}}\right)\text{Myr}.$$(78)

Clearly, the YORP timescale is shorter than τ_{imp,ire–ori}, although its specific value exhibits considerable variation across different asteroids. Consequently, it is widely accepted that the YORP effect primarily governs the rotational evolution of small objects, while collisions play a dominant role in larger objects. The classic static YORP model – which assumes an invariable YORP torque until the asteroid spins up to disruption or spins down to a quasi-static rotation – has been used to study the long-term rotational evolution of asteroids (Rubincam 2000; Pravec et al. 2005; Bottke et al. 2015). A more intricate model, referred to as the ‘stochastic YORP model’ (Bottke et al. 2015), takes into account the resetting of the YORP torque caused by collisions, which arises from the high sensitivity of the YORP effect to surface morphology. Although a suggested timescale of 1 Myr has been proposed for YORP reset (Bottke et al. 2015), a quantitative assessment of the torque changes resulting from craters is yet to be explored. The CYORP effect offers a powerful tool for investigating the stochastic YORP effect. While a comprehensive examination of the stochastic YORP effect is beyond the scope of this paper, we present an example of integrating the CYORP effect into the static YORP model.

In our simulation, a random YORP coefficient is assigned within the range of −0.005 to 0.005, with the coefficient’s sign matching that of the torque. The CYORP torque is introduced specifically when a collision event takes place. The timescale for the impact by an asteroid with the size R_imp is $${\tau }_{\text{imp}}=\frac{1}{P_{\text{i}}\pi R_{\text{Antonia }}^2\Delta N\left(R>R_{\text{imp}}\right)},$$(79)

where $$N\left(R>R_{\text{imp}}\right)=C_R{\left(\frac{R_{\text{imp}}}{1 \text{km}}\right)}^{-b_R}.$$(80)

Here, P_i = 2.85 × 10⁻¹⁸ km⁻² yr⁻¹ is the intrinsic collision probability, C_R = 6 × 10⁵, and b_R = 2.2 (Holsapple 2022). In the strength regime formulation, the crater produced by an impactor with the size of R_imp has a size of $$R_{\text{crater }}=1.306R_{\text{imp }}\left({\frac{{\rho }_{\text{imp }}}{{\rho }_{\text{ast }}}}^{0.4}\right){\left(\frac{Y}{\rho v_{\text{imp }}}\right)}^{-0.205},$$(81)

with Y = 100 Pa and υ_imp = 5.3 km s⁻¹ .In the gravity regime, $$R_{\text{crater }}=0.59{\left(\frac{{\rho }_{\text{ast }}}{m_{\text{imp }}}\right)}^{-1/3}{\left(\frac{{\rho }_{\text{ast }}}{{\rho }_{\text{imp }}}\right)}^{0.00138}{\left(\frac{g_{\text{ast }}{R}_{\text{ast }}}{v_{\text{imp }}}\right)}^{-0.17},$$(82)

where 𝑔_ast is the surface gravity of the asteroid. In each time step (~10³ yr), we calculate the numbers of impact craters of different sizes, according to Eqs. (79) and (81). We then assign each crater with a random surface element of the polyhedron model of the asteroid Antonia, after which we can calculate the CYORP torque. Finally, we add the CYORP torque to the normal YORP torque directly calculated from the shape model in order to obtain the total YORP torque. The spin rate evolves following $$\dot{\omega }=\frac{T_z}{I},$$(83)

with I being the maximum moment of inertia and T_z being the torque component that is along a spin vector. The obliquity evolves according to $$\dot{ϵ}=\frac{T_{ϵ}}{I\omega },$$(84)

where T_ϵ is the torque component that changes the obliquity.

There exist two possible end states in a YORP cycle: either the asteroid’s rotation slows down until it reaches a quasi-non-rotational state, or it accelerates to the spin threshold for shape change or disruption with a period of approximately 2.2 h.

Upon completing a YORP cycle, the asteroid’s rotational state is updated by assigning a new random rotational speed and obliquity. The impact of introducing CYORP torques can be observed in the evolution of a 10 km asteroid, as depicted in Fig. 11. Notably, significant differences arise when considering the inclusion of CYORP torques.

Nonetheless, the rotational evolution of asteroids currently lacks a standardised model. Some models propose that after spinning down to a non-rotational state, the asteroid’s spin rate is assigned a new random value within a specified range (Hanus et al. 2011; Bottke et al. 2015), while some assume it continues to spin up under the YORP effect (Pravec et al. 2008; Marzari et al. 2020). By selecting an initial spin rate for a new rotational state, Holsapple (2022) reproduces the spin evolution without the YORP effect. Hence, rather than attempting to address the entire complexity of the problem, our objective in this study is to present an illustrative example of the interaction between CYORP and the conventional YORP effect. Furthermore, we underscore the significance of the CYORP effect in the long-term rotational evolution of asteroids. A comprehensive investigation of the rotational evolution of asteroid groups is left for future research.

Fig. 10 YORP torque damped by the CYORP effect in the case of asteroid (272) Antonia. The spin component and the obliquity component are shown in the left and middle panels, respectively. The right panel shows the CYORP coefficient distribution over the asteroid’s surface.

Fig. 11 Evolution of the spin rate (left panel) and the obliquity (right panel) of a 10 km synthetic asteroid. In the presence of the static YORP torque (blue line), the asteroid gradually decelerates until it reaches a quasi-non-rotational state. Subsequently, we impose a new rotational state on the asteroid by assigning random values of rotational speed and obliquity. Conversely, when incorporating the CYORP torque (red line), the asteroid follows a different path, exhibiting random fluctuations in its spin rate due to the occurrence of impacts, creating new craters that lead to changes in the CYORP torque. As a result, the asteroid experiences intermittent transitions between spin up and spin down.

7 Summary and conclusions

The YORP effect is a thermal torque produced by radiation from the irregular surface of the asteroid. It has been demonstrated that this effect is highly sensitive to surface topology (Statler 2009; Breiter et al. 2009), including small-scale roughness (Rozitis & Green 2012), boulders (Golubov & Krugly 2012), and craters (Zhou et al. 2022). In this study, we developed a semi-analytical model for calculating the temperature field of a crater, which accounts for the effects of self-sheltering, self-radiation, self-scattering, and non-zero thermal conductivity. Using this model, we investigated the crater-induced YORP (CYORP) effect in a computationally efficient manner (about three orders of magnitude faster than the numerical method), allowing for a comprehensive exploration of the functional dependence of the CYORP effect and its incorporation into the rotational and orbital evolution of asteroids. The main results and conclusions of this study can be summarised as follows.

Our semi-analytical model for the CYORP effect is valid in the high-thermal-conductivity regime (K > 0.1 Wm⁻¹ K⁻¹). This suggests that the model is suitable for application to materials such as solid basalt and metal, which are usually beneath the regolith on asteroid surfaces but may be exposed to sunlight due to the formation of deep craters.

The CYORP effect is significant when the crater depth-to-diameter ratio is greater than 0.05. The self-modification effects of a concave structure, including the self-sheltering effect, self-radiation effect, and self-scattering effect, are stronger with a higher depth-to-diameter ratio. For concave structures with a depth-to-diameter ratio of smaller than 0.05, the surface can be treated as a convex shape without introducing significant inaccuracies. The typical value of the CYORP coefficient for the spin component is 0.01, which is insensitive to the thermal parameter, while the obliquity component decreases from 0.004 for basalt to 0.001 for metal. For a z-axis symmetric shape (e.g. a spinning top shape), the spin component of the CYORP torque vanishes while the obliquity component survives, which implies that the spin acceleration of such symmetric shapes does not change significantly under the effect of crater formation.

Using our semi-analytical method, we confirm that the YORP torque can be damped by the surface roughness, which was first discovered by Rozitis & Green (2012). The fast computation of our semi-analytical model allows us to consider more flexible configurations of surface roughness, such as a space-varying roughness distribution, roughness on components of binary asteroids, and so on.

The magnitude and directional change of the YORP coefficient at each impact are assessed for the first time using our CYORP model. While a complete investigation of the spin evolution of asteroids is left for future work, we show that rotational evolution can be severely affected by collisions.

Acknowledgements

We acknowledge support from the Université Côte d’Azur. We thank Yun Zhang, Xiaoran Yan, and Marco Delbo for the useful discussion. Wen-Han Zhou would like to acknowledge the funding support from the Chinese Scholarship Council (No. 202110320014). Patrick Michel acknowledges funding support from the French space agency CNES and from the European Union’s Horizon 2020 research and innovation program under grant agreement no. 870377 (project NEO-MAPP).

References

All Tables

All Figures

	Fig. 1 Three coordinate systems in this paper: coordinate system oxyz for calculating the illuminated domain in the crater, PABC for calculating the effective radiative force of an arbitrary surface element, and OXYZ for averaging the YORP torque over the spin and orbital motion.
In the text

Fig. 2 Radiative force of the total crater averaged over a rotational period (8 h by default) obtained from the analytical method (solid lines) and the numerical method (dashed lines), as a function of the crater colatitude α for K = 2.65 Wm−1 K−1 (left panel) and K = 0.1 Wm−1 K−1 (right panel). The x, y, and 𝓏 components of the radiative force are shown in blue, red, and green, respectively. The distance of the crater here is 1 au from the Sun.

In the text

	Fig. 3 Same as Fig. 2 but for K = 0 Wm−1 K−1 in the analytical model and K = 0.001 Wm−1 K−1 in the numerical model.
In the text

	Fig. 4 Ratio of the radiative force including self-modification effects to that without self-modification effects, as a function of crater depth-to-diameter ratio and accounting for different self-modification effects.
In the text

	Fig. 5 Spin component (left panel) and the obliquity component (right panel) of the CYORP coefficient as a function of the asteroid obliquity e, accounting for different depth-to-diameter ratios, which are denoted by different colours.
In the text

	Fig. 6 Spin component (left panel) and the obliquity component (right panel) of the CYORP coefficient as a function of the asteroid obliquity e and the crater colatitude α(> ∆), for K = 2.65Wm−1 K−1. Here, d/D0 = 0.16 and δ = ∆ = π/4.
In the text

	Fig. 7 Same as Fig. 6 but for K = 40 Wm−1 K−1.
In the text

	Fig. 8 Spin component (left panel) and the obliquity component (right panel) of the CYORP coefficient as a function of ∆(<α) and δ, for K = 2.65 Wm−1 K−1. Here, d/D0 = 0.16 and a = ϵ = π/6.
In the text

	Fig. 9 Obliquity component (right panel) of the CYORP coefficient as a function of the asteroid obliquity e and the crater colatitude α for K = 2.65 Wm−1 K−1. Here, d/D0 = 0.16 and δ = ∆ = 0.
In the text

	Fig. 10 YORP torque damped by the CYORP effect in the case of asteroid (272) Antonia. The spin component and the obliquity component are shown in the left and middle panels, respectively. The right panel shows the CYORP coefficient distribution over the asteroid’s surface.
In the text

Fig. 11 Evolution of the spin rate (left panel) and the obliquity (right panel) of a 10 km synthetic asteroid. In the presence of the static YORP torque (blue line), the asteroid gradually decelerates until it reaches a quasi-non-rotational state. Subsequently, we impose a new rotational state on the asteroid by assigning random values of rotational speed and obliquity. Conversely, when incorporating the CYORP torque (red line), the asteroid follows a different path, exhibiting random fluctuations in its spin rate due to the occurrence of impacts, creating new craters that lead to changes in the CYORP torque. As a result, the asteroid experiences intermittent transitions between spin up and spin down.

In the text