# A semi-analytical thermal model for craters with application to the crater-induced YORP effect

Wen-Han Zhou<fnref n="1" />, and Patrick Michel<fnref n="1" />,<fnref n="2" />

Received 23 May 2023 / Accepted 4 December 2023

## ABSTRACT

**Context.** The YORP effect is the thermal torque generated by radiation from the surface of an asteroid. The effect is sensitive to surface topology, including small-scale roughness, boulders, and craters.

**Aims.** The aim of this paper is to develop a computationally efficient semi-analytical model for the crater-induced YORP (CYORP) effect that can be used to investigate the functional dependence of this effect.

**Methods.** This study linearizes the thermal radiation term as a function of the temperature in the boundary condition of the heat conductivity, and obtains the temperature field in a crater over a rotational period in the form of a Fourier series, accounting for the effects of self-sheltering, self-radiation, and self-scattering. By comparison with a numerical model, we find that this semi-analytical model for the CYORP effect works well for $K > 0.1 \text{ W m}^{-1} \text{ K}^{-1}$. This semi-analytical model is computationally three-orders-of-magnitude more efficient than the numerical approach.

**Results.** We obtain the temperature field of a crater, accounting for the thermal inertia, crater shape, and crater location. We then find that the CYORP effect is negligible when the depth-to-diameter ratio is smaller than 0.05. In this case, it is reasonable to assume a convex shape for YORP calculations. Varying the thermal conductivity yields a consistent value of approximately 0.01 for the spin component of the CYORP coefficient, while the obliquity component is inversely related to thermal inertia, declining from 0.004 in basalt to 0.001 in metal. The CYORP spin component peaks at an obliquity of 0°, 90°, or 180°, while the obliquity component peaks at an obliquity of around 45° or 135°. For a z-axis symmetric shape, the CYORP spin component vanishes, while the obliquity component persists. Our model confirms that the total YORP torque is damped by a few tens of percent by uniformly distributed small-scale surface roughness. Furthermore, for the first time, we calculate the change in the YORP torque at each impact on the surface of an asteroid explicitly and compute the resulting stochastic spin evolution more precisely.

**Conclusions.** This study shows that the CYORP effect due to small-scale surface roughness and impact craters is significant during the history of asteroids. The semi-analytical method that we developed, which benefits from fast computation, offers new perspectives for future investigations of the YORP modeling of real asteroids and for the complete rotational and orbital evolution of asteroids accounting for collisions. Future research employing our CYORP model may explore the implications of space-varying roughness distribution, roughness in binary systems, and the development of a comprehensive rotational evolution model for asteroid groups.

**Key words.** minor planets, asteroids: general

## 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; Vokrouhlický 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 (Rožek 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 (Cuk & 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

<footnote n="1">¹ Université Côte d'Azur, Observatoire de la Côte d'Azur, CNRS, Laboratoire Lagrange, 96 Bd de l'Observatoire, 06304 Nice, France
e-mail: wenhan.zhou@oca.eu</footnote>

<footnote n="2">² University of Tokyo, Department of Systems Innovation, School of Engineering, Tokyo, Japan</footnote>
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}}. \qquad (1)$$

In general, it takes the form of the following scaling rule with the radius of the crater $R_0$ and of the asteroid radius $R_{\text{ast}}$:

$$T_{\text{CYORP}} = W \frac{\Phi}{c} R_0^2 R_{\text{ast}}, \qquad (2)$$

where $\Phi$ 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}, \qquad (3)$$

with two boundary conditions,

$$K \frac{\partial u}{\partial z} \Big|_{z=0} = F(t) - e \sigma u^4 \Big|_{z=0}, \qquad (4)$$

$$K \frac{\partial u}{\partial z} \Big|_{z \to \infty} = 0, \qquad (5)$$

and a periodic condition,

$$u|_{t=2\pi/\omega} = u|_{t=0}, \qquad (6)$$

where $t$ is the time, $z$ is the depth of the crater, $\omega$ is the angular velocity, $K$ is the thermal conductivity, $C$ is the specific heat capacity, $\rho$ is the bulk density of the asteroid, $e$ is the emissivity, and $\sigma$ is the Stefan-Boltzmann constant. In the following, we use the spin angle $\beta = \omega t$ to replace time for simplicity. The effect of the seasonal wave is marginal and is ignored in this work.
<figure><img src="image_1.png" /><figcaption>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.</figcaption></figure>

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 $\theta$ and the azimuthal angle $\phi$, the solution of the temperature is the real part of the expression

$$u(\theta, \phi, \beta) = u_0(\theta, \phi) + \sum_{n=1}^{\infty} u_n(\theta, \phi)e^{in\beta}, \quad (7)$$

with

$$u_0(\theta, \phi) = \left( \frac{F_0(\theta, \phi)}{e\sigma} \right)^{1/4}, \quad (8)$$

$$u_n(\theta, \phi) = \frac{F_n(\theta, \phi)}{4e\sigma u_0(\theta, \phi)^3} \frac{e^{iJ_n}}{\sqrt{2\Theta_n^2 + 2\Theta_n + 1}}, \quad (9)$$

where

$$l_n = \sqrt{\frac{n\omega}{2M}}, \quad (10)$$

$$\Theta_n = \frac{Kl_n}{4e\sigma u_0^3}, \quad (11)$$

$$\tan J_n = -\frac{\Theta_n}{\Theta_n + 1}, \quad (12)$$

with $M = K/\rho 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, \phi) = \sum_{n=0}^{\infty} F_n(\theta, \phi)e^{in\beta}. \quad (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(\theta, \phi, \beta)$, radiation from the crater itself $H(\theta, \phi, \beta)$, and the scattering flux from the crater $G(\theta, \phi, \beta)$, which leads us to

$$\begin{aligned} F(\theta, \phi) &= (1 - A)(E(\theta, \phi) + H(\theta, \phi) + G(\theta, \phi)) \\ &= (1 - A) \sum_{n=0}^{\infty} (E_n(\theta, \phi) + H_n(\theta, \phi) + G_n(\theta, \phi))e^{in\beta}, \end{aligned} \quad (14)$$

where $E_n(\theta, \phi, \beta)$ and $H_n(\theta, \phi, \beta)$ denote the nth-order Fourier modes of $E(\theta, \phi, t)$ and $H(\theta, \phi, 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, \phi) = \left( \frac{(1 - A)(E_0(\theta, \phi) + H_0(\theta, \phi) + G_0(\theta, \phi))}{e\sigma} \right)^{1/4}, \quad (15)$$

$$u_n(\theta, \phi) = \frac{(1 - A)(E_n(\theta, \phi) + H_n(\theta, \phi) + G_n(\theta, \phi))}{4e\sigma u_0(\theta, \phi)^3} \frac{e^{iJ_n}}{\sqrt{2\Theta_n^2 + 2\Theta_n + 1}}. \quad (16)$$

Therefore, to solve the temperature $u(\theta, \phi, \beta)$, we need to obtain the Fourier series of $E(\theta, \phi, \beta)$, $H(\theta, \phi, \beta)$, and $G(\theta, \phi, \beta)$.

### 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_1$ 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 $\mathbf{r}_0$ (denoted by $\alpha$ and $\beta$), the crater normal vector $\mathbf{n}_0$ (described by $\alpha, \beta, \delta$, and $\Delta$), and the solar position $\mathbf{s}$ (described by $\epsilon$ and $\eta$), as shown in the right panel of Fig. 1. As the polar angle is between 0° and 180° by definition, we limit $0° < \alpha < 180°$ and $0° < \alpha - \Delta < 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 $\mathbf{n}_0$, which is also the symmetry axis of the spherical crater. The direction of axis $oy$ is along $\mathbf{e}_{OZ} \times \mathbf{e}_{OZ}$ and $ox$ follows the right-hand rule. In this coordinate system, a crater with a depth of $h$ can be defined as

$$Z := \{(x, y, z) \in \mathbb{R}^3 | x^2 + y^2 + z^2 = R_1, z \ge R_1 \cos \theta_0\}, \quad (17)$$

$$\text{with } \cos \theta_0 = (R_1 - h)/R_1.$$

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 $\mathbf{e}_{OZ} \times \mathbf{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 $\mathbf{r}_P = (\sin \theta \cos \phi, \sin \theta \sin \phi, \cos \theta)$. The unit position vector of the Sun in the coordinate system OXYZ is described by

$$s_{OXYZ} = (\cos \eta, \cos \epsilon \sin \eta, \sin \epsilon \sin \eta), \quad (18)$$

where $\epsilon$ is the obliquity and $\eta$ is the angle of orbital motion. The transform matrix between the coordinate systems $oxyz$ and OXYZ is set to

$$\mathcal{R} = \begin{pmatrix} \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{pmatrix}. \quad (19)$$

Therefore, the coordinates of the vector $\mathbf{s}$ in the coordinate system $oxyz$ is

$$s_{oxyz} = \begin{pmatrix} -\sin \alpha \sin \epsilon \sin \eta + \cos \alpha \cos \beta \cos \eta + \cos \alpha \sin \beta \cos \epsilon \sin \eta \\ \sin \beta \cos \eta - \cos \beta \cos \epsilon \sin \eta \\ -\sin \alpha \cos \beta \cos \eta - \sin \alpha \sin \beta \cos \epsilon \sin \eta - \cos \alpha \sin \epsilon \sin \eta \end{pmatrix}. \quad (20)$$

On the other hand, we can use the angle $\lambda$ and $\phi'$ to represent $s_{oxyz}$:

$$s_{oxyz} = (\sin \lambda \cos \phi', \sin \lambda \sin \phi', -\cos \lambda), \quad (21)$$

such that

$$\cos \lambda = \sin \alpha \cos \beta \cos \eta + \sin \alpha \sin \beta \cos \epsilon \sin \eta + \cos \alpha \sin \epsilon \sin \eta, \quad (22)$$

and

$$\tan \phi' = \frac{\sin \alpha \sin \epsilon \sin \eta - \cos \alpha \cos \beta \cos \eta - \cos \alpha \sin \beta \cos \epsilon \sin \eta}{\sin \beta \cos \eta + \cos \beta \cos \epsilon \sin \eta}. \quad (23)$$

The absorbed radiation flux is

$$\begin{aligned} E(\theta, \phi) &= (1 - A)\Phi H(\cos \lambda)H(w)(-\mathbf{r}_P \cdot \mathbf{s}) \\ &= (1 - A)\Phi H(\cos \lambda)H(w) \cdot [m_1 \cos(\beta - \beta_1) + m_2] \end{aligned} \quad (24)$$

where

$$w = \cos 2\lambda \cos \theta + \sin \theta_0 - \sin 2\lambda \sin \theta \cos(\phi - \phi'), \quad (25)$$

$$\tan \beta_1 = \frac{(\sin \alpha \cos \theta - \cos \alpha \sin \theta \cos \phi) \cos \epsilon \sin \eta - \sin \theta \sin \phi \cos \eta}{(\sin \alpha \cos \theta - \cos \alpha \sin \theta \cos \phi) \cos \eta + \cos \epsilon \sin \eta \sin \theta \sin \phi}, \quad (26)$$

$$\begin{aligned} m_1 = (\cos^2 \eta + \cos^2 \epsilon \sin^2 \eta)^{1/2} & \cdot (\cos^2 \theta \sin^2 \alpha - \frac{\cos \phi \sin 2\alpha \sin 2\theta}{2}) \\ & + (\cos^2 \alpha \cos^2 \phi + \sin^2 \phi) \sin^2 \theta)^{1/2}, \end{aligned} \quad (27)$$

$$m_2 = \cos \alpha \cos \theta \sin \epsilon \sin \eta + \cos \phi \sin \alpha \sin \epsilon \sin \eta \sin \theta. \quad (28)$$

Here, $H$ is the Heaviside function defined by

$$H(x) := \begin{cases} 1, & x > 0 \\ 0, & \text{else.} \end{cases} \quad (29)$$

Using the substitution $\beta' = \beta - \beta_1$, the absorbed radiation flux has the form

$$E(\theta, \phi) = \begin{cases} (1 - A)\Phi(m_1 \cos \beta' + m_2), & \beta'_{\min} < \beta < \beta'_{\max} \\ 0, & \text{else.} \end{cases} \quad (30)$$

Expanding Eq. (30) in Fourier series, we obtain

$$E(\theta, \phi) = E_0 + \frac{(1 - A)\Phi}{\pi} \sum_{n=1}^{\infty} [C_n \cos(n\beta') + S_n \sin(n\beta')], \quad (31)$$

with

$$E_0 = \frac{(1 - A)\Phi}{2\pi} (m_1 \sin \beta + m_2\beta) |\beta'_{\min}|^{\beta'_{\max}}, \quad (32)$$

$$C_n = \int_{\beta'}^{\beta'_{\min}} (m_1 \cos \beta'_{\max} + m_2) \cos(n\beta') d\beta' = \begin{cases} \left( \frac{m_1 \sin 2\beta' + 2m_1\beta'}{4} + m_2 \sin \beta \right) |\beta'_{\max}|, & n = 1 \\ \left( \frac{m_1 n \cos \beta' \sin n\beta' - m_1 \cos n\beta' \sin \beta'}{n^2 - 1} + \frac{m_2 \sin(n\beta')}{n} \right) |\beta'_{\max}|, & n > 1, \end{cases} \quad (33)$$
and

$$
\begin{align}
S_n &= \int_{\beta'}^{\beta'_{\max}} (m_1 \cos \beta'_{\max} + m_2) \sin(n\beta') d\beta' \nonumber \\
&= \begin{cases} -\left(\frac{m_1 \cos \beta'^2}{2} + m_2 \cos \beta'\right) |\beta'_{\max}|, & n=1 \\ \left(\frac{m_1 n \cos \beta' \cos n\beta' + m_1 \sin \beta' \sin n\beta'}{1-n^2} - \frac{m_2 \cos n\beta'}{n}\right) |\beta'_{\max}|, & n>1. \end{cases} \tag{34}
\end{align}
$$

The nth coefficient of the Fourier series of E(θ, φ) is

$$
E_n = \sqrt{C_n^2 + S_n^2} \cdot e^{i\Phi_n}, \quad (35)
$$

with

$$
\tan \Phi_n = -\frac{S_n}{C_n}. \tag{36}
$$

2.1.4. Self-heating effect

Due to the concave structure of the crater, the surface element in
the crater also receives photons emitted by the crater itself, 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

$$
\mathbf{r}_1 = R_1(\sin \theta_1 \cos \phi_1, \sin \theta_1 \sin \phi_1, \cos \theta_1), \quad (37)
$$

$$
\mathbf{r}_2 = R_1(\sin \theta_2 \cos \phi_2, \sin \theta_2 \sin \phi_2, \cos \theta_2). \quad (38)
$$

The displacement from dS₁ to dS₂ is

$$
\mathbf{r}_{1,2} = \mathbf{r}_2 - \mathbf{r}_1. \tag{39}
$$

The incident angle $\zeta_1$ and the emission angle $\zeta_2$ are defined as

$$
\cos \zeta_1 = -\mathbf{r}_1 \cdot \mathbf{r}_{1,2} / R_1 r_{1,2}, \quad (40)
$$

$$
\cos \zeta_2 = \mathbf{r}_2 \cdot \mathbf{r}_{1,2} / R_1 r_{1,2}, \quad (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 \cos \zeta_1 \cos \zeta_2}{\pi \mathbf{r}_{1,2}^2} dS_2. \quad (42)
$$

Substituting Eqs. (37)–(41) into Eq. (42), we obtain

$$
H_{1,2} = \frac{e \sigma u_2^4}{4 \pi R_1^2} dS_2. \quad (43)
$$

For an arbitrary point, the radiation flux caused by the whole
crater is

$$
H(\theta, \phi) = \int_Z \frac{e \sigma u(\theta', \phi')^4}{4\pi} \sin \theta' d\theta' d\phi', \quad (44)
$$

where $\mathcal{Z}$ is the crater surface, and is defined as

$$
\mathcal{Z} := \{(x, y, z) \in \mathbb{R}^3 | r = R_1, \theta \in (0, \pi/2 - \theta_0), \phi \in (0, 2\pi)\}. \quad (45)
$$

Similarly, we obtain the self-scattering term:

$$
G(\theta, \phi) = \int_Z \frac{AE(\theta', \phi')}{4\pi} \sin \theta' d\theta' d\phi'. \quad (46)
$$

Therefore, both $H(\theta, \phi)$ and $G(\theta, \phi)$ can be expressed in terms of $u(\theta, \phi)$ and $E(\theta, \phi)$. As $E(\theta, \phi)$ is derived in Sect. 2.1.3, the only unknown variable is the temperature $u(\theta, \phi)$, 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(\theta, \phi)$. 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{align}
\sigma u_0^4(\theta, \phi) &= (1-A)(E_0(\theta, \phi) + G_0(\theta, \phi) + H_0(\theta, \phi)) \nonumber \\
\sigma u_n(\theta, \phi) u_0^3(\theta, \phi) \frac{2\sqrt{2\Theta_n^2 + 2\Theta_n + 1}}{(1-A)e^{iJ_n}} &= E_n(\theta, \phi) + G_n(\theta, \phi) + H_n(\theta, \phi), \tag{47}
\end{align}
$$

with

$$
\begin{align}
G_n(\theta, \phi) &= \int_Z \frac{AE_n(\theta', \phi')}{4\pi} \sin\theta' d\theta' d\phi', \nonumber \\
H_0(\theta, \phi) &= \int_Z \frac{e\sigma u_0(\theta', \phi')^4}{4\pi} \sin\theta' d\theta' d\phi', \tag{48} \\
H_n(\theta, \phi) &= \int_Z \frac{e\sigma u_0(\theta', \phi')^3 u_n(\theta', \phi')}{\pi} \sin\theta' d\theta' d\phi'. \nonumber
\end{align}
$$

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, \phi) = g(\theta, \phi) + C \int_Z f(\theta', \phi') \sin \theta' d\theta' d\phi'. \quad (49)
$$

By setting

$$
k = \int_Z f(\theta', \phi') \sin \theta' d\theta' d\phi', \quad (50)
$$

we have

$$
f(\theta, \phi) = g(\theta, \phi) + C \cdot k. \quad (51)
$$

Substituting Eq. (51) into Eq. (50), we obtain

$$
k = \int_Z (g(\theta', \phi') + C \cdot k) \sin \theta' d\theta' d\phi'. \quad (52)
$$

Rearranging Eq. (52), we find $k$ is

$$
k = \frac{\int_Z g(\theta', \phi') \sin \theta' d\theta' d\phi'}{1 - C \int_Z \sin \theta' d\theta' d\phi'}, \quad (53)
$$

with which $f(\theta, \phi)$ is solved out by Eq. (51).

2.2.2. Solution for temperature

In the case of $u_0$,

$$
f(\theta, \phi) = e \sigma u_0^4 / (1 - A), \quad (54)
$$

$$
g(\theta, \phi) = E_0(\theta, \phi) + G_0(\theta, \phi), \quad (55)
$$

$$
C = \frac{1}{4\pi}, \tag{56}
$$
and in the case of $u_n$,

$$f(\theta, \phi) = e\sigma u_0^3 u_n \frac{2\sqrt{2\Theta_n^2 + 2\Theta_n + 1}}{(1 - A)e^{iJ_n}}, \quad (57)$$

$$g(\theta, \phi) = E_n(\theta, \phi) + G_n(\theta, \phi), \quad (58)$$

$$C = \frac{e^{iJ_n}}{2\pi \sqrt{2\Theta_n^2 + 2\Theta_n + 1}}. \quad (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}}. \quad (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}} = \mathbf{r}_{0,\text{OXYZ}} \times f_{\text{OXYZ}}. \quad (61)$$

Here, $f_{\text{OXYZ}}$ is the radiative force, and

$$f_{\text{OXYZ}} = \int_W \frac{e\sigma T(\theta, \phi)^4}{c} \mathbf{n}_f(\theta, \phi) dS, \quad (62)$$

where $\mathbf{n}_f$ is the corrected force direction vector for each surface element. If there is no shelter around the surface element, $\mathbf{n}_f$ is equal to the surface unit normal vector $\mathbf{n}_0$. 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(\theta, \phi)$, the radiative force is

$$f = - \int_{H} \frac{\Phi}{\pi c} \cos \theta' \begin{pmatrix} \sin \theta' \cos \phi' \\ \sin \theta' \sin \phi' \\ \cos \theta' \end{pmatrix} \sin \theta' d\theta' d\phi'. \quad (63)$$

Without sheltering (e.g. for convex asteroids), $H$ is replaced by the hemisphere (i.e. $\theta \in (0, \pi/2)$, $\phi \in (0, 2\pi)$). In this case, the force is reduced to $2\Phi\mathbf{n}_0/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 $\mathbf{n}_0$ is

$$E_{\text{ground}} = (1 - A)\Phi H(\cos \lambda) \cos \lambda \\ = \begin{cases} (1 - A)\Phi(m_1 \cos \beta' + m_2), & \beta'_{\min} < \beta' < \beta'_{\max} \\ 0, & \text{else.} \end{cases} \quad (64)$$

Here, $\beta' = \beta - \beta_1$, with

$$\tan \beta_1 = \frac{\sin \alpha \cos \epsilon \sin \eta}{\sin \alpha \cos \eta}, \quad (65)$$

$$m_1 = (\sin^2 \alpha (\cos^2 \eta + \cos^2 \epsilon \sin^2 \eta))^{1/2}, \quad (66)$$

$$m_2 = \cos \alpha \sin \epsilon \sin \eta, \quad (67)$$

and $\beta_{\min}$ and $\beta_{\max}$ are the negative and positive values of $\arccos(-m_2/m_1)$, 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 \sim u_0^3 \sum u_n e^{\beta n}$). This assumption is equivalent to a high value of the thermal parameter, which is defined as

$$\Gamma = \frac{\sqrt{C \rho \omega \kappa}}{(\epsilon \sigma)^{1/4} (1 - A)^{3/4} \Phi^{3/4}}. \quad (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:

$$c_n u_{i,j+1,k+1} - (2c_n + 1)u_{i,j,k+1} + c_n u_{i,j-1,k+1} = -c_n u_{i,j+1,k} + (2c_n - 1)u_{i,j,k} - c_n u_{i,j-1,k}. \quad (69)$$

Here, $u_{i,j,k}$ represents the temperature at the depth of $(j-1)\delta z$ below the $i$th facet at the $k$th 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_n = \frac{a\delta t}{2(\delta x)^2}. \quad (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)(E_{i,k} + H_{i,k} + S_{i,k}) - \sigma u_{i,1,k+1}^4 = K \frac{u_{i,1,k+1} - u_{i,2,k}}{\delta z}, \quad (71)$$

which can be solved by a Newtonian-Raphson iterative method. The solar flux $E_{i,k}$ on the $i$th facet at the $k$th time step is

$$E_{i,k} = \begin{cases} \Phi_0 \mathbf{n}_i \cdot \mathbf{s}_k H(\mathbf{n}_i \cdot \mathbf{s}_k), & \text{unsheltered} \\ 0, & \text{sheltered}, \end{cases} \quad (72)$$

where $\Phi_0$ is 1364 W m$^{-2}$ at the distance of 1 au to the Sun. Here, $\mathbf{n}_i$ is the normal vector of the $i$th facet and $\mathbf{s}_k$ is the unit solar position vector at the $k$th time step. Whether or not the facet is sheltered is judged according to the projections of other facets on the plane of the $i$th facet along the solar position vector $\mathbf{s}_k$. The self-radiation term $H_{i,k}$ is the sum of radiation from other facets:

$$H_{i,k} = \sum_{i' \neq i} e\sigma u_{i',1,k}^4 \frac{-(\mathbf{n}_i \cdot \mathbf{r}_{i,i'}) (\mathbf{n}_{i'} \cdot \mathbf{r}_{i,i'})}{\pi r_{i,i'}^2} S_{i'}, \quad (73)$$

where $\mathbf{r}_{i,i'}$ is the vector from the centre of the $i$th facet to the centre of $i'$th facet, and $S_{i'}$ is the area of the $i$th facet. The scattering term $S_{i,k}$ is given by

$$S_{i,k} = \sum_{i' \neq i} A E_{i'} \frac{-(\mathbf{n}_i \cdot \mathbf{r}_{i,i'})(\mathbf{n}_{i'} \cdot \mathbf{r}_{i,i'})}{\pi r_{i,i'}^2} S_{i'}. \quad (74)$$
<figure><img src="image_2.png" /><figcaption>**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 $\alpha$ for $K = 2.65 \text{ W m}^{-1} \text{ K}^{-1}$ (left panel) and $K = 0.1 \text{ W m}^{-1} \text{ K}^{-1}$ (right panel). The $x$, $y$, and $z$ 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.</figcaption></figure>

The second boundary condition (Eq. (5)) translates into

$$u_{i,j_{\max},k+1} = u_{i,j_{\max}-1,k+1}. \qquad (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 \text{ 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}\delta t \sim 20$ spin periods. We divided the considered crater into about 1000 facets and solve the temperature for each facet using the above method.

## 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 \text{ W mK}^{-1}$ to $1 \text{ W mK}^{-1}$, while the thermal conductivity of metallic asteroids is generally much higher, in the range of $20–50 \text{ W mK}^{-1}$. 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 $\epsilon = 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

<table><caption>**Table 1.** Thermal parameter for three typical materials on asteroids, taken from Farinella et al. (1998).</caption><thead><tr><th></th><th>ρ (kg m<sup>-3</sup>)</th><th>K (W m<sup>-1</sup> K<sup>-1</sup>)</th><th>C (J kg<sup>-1</sup> K<sup>-1</sup>)</th></tr></thead><tbody><tr><td>Regolith</td><td>1500</td><td>0.0015</td><td>680</td></tr><tr><td>Solid basalt</td><td>3500</td><td>2.65</td><td>680</td></tr><tr><td>Solid iron</td><td>8000</td><td>40</td><td>500</td></tr></tbody></table>

decreases to $0.1 \text{ W m}^{-1} \text{ K}^{-1}$. 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 \text{ W m}^{-1} \text{ K}^{-1}$, 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
<figure><img src="image_3.png" /><figcaption>Fig. 3. Same as Fig. 2 but for K = 0 W m⁻<fnref n="1" /> K⁻<fnref n="1" /> in the analytical model and K = 0.001 W m⁻<fnref n="1" /> K⁻<fnref n="1" /> in the numerical model.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>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.</figcaption></figure>

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 < 0.05$, the force increase is less than 1%. Therefore, for those shallow concave structures with $h/D_0 < 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.

## 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$ W m⁻<fnref n="1" /> K⁻<fnref n="1" />, and $α = π/2$ are used. Higher depth-to-diameter ratios correspond to larger spin components of the CYORP torque. A crater with $h/D_0 < 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_z$ 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 asteroids<fnref n="1" />. The ratio of the CYORP torque to the normal YORP torque scales as

$$ \frac{T_{CYORP}}{T_{NYORP}} = \frac{W_{CYORP}S_{crater}}{W_{CYORP}S_{asteroid}} \qquad (76) $$

<footnote n="1">¹ See the definition by Vokrouhlický & Čapek (2002).</footnote>
<figure><img src="image_5.png" /><figcaption>**Fig. 5.** Spin component (left panel) and the obliquity component (right panel) of the CYORP coefficient as a function of the asteroid obliquity $\epsilon$, accounting for different depth-to-diameter ratios, which are denoted by different colours.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Fig. 6.** Spin component (left panel) and the obliquity component (right panel) of the CYORP coefficient as a function of the asteroid obliquity $\epsilon$ and the crater colatitude $\alpha(> \Delta)$, for $K = 2.65\,\text{W m}^{-1} \text{K}^{-1}$. Here, $d/D_0 = 0.16$ and $\delta = \Delta = \pi/4$.</figcaption></figure>

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.

### 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 \, \text{W m}^{-1} \text{K}^{-1}$ 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 $\delta$ and $\Delta$

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

Zhou et al. (2022) demonstrates that the CYORP torque vanishes for $\delta = 0$. However, in the presence of finite thermal inertia, the obliquity component of the CYORP torque arises while
<figure><img src="image_7.png" /><figcaption>Fig. 7. Same as Fig. 6 but for $K = 40 \text{ W m}^{-1} \text{ K}^{-1}$.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 8. Spin component (left panel) and the obliquity component (right panel) of the CYORP coefficient as a function of $\Delta(<\alpha)$ and $\delta$, for $K = 2.65 \text{ W m}^{-1} \text{ K}^{-1}$. Here, $d/D_0 = 0.16$ and $\alpha = \epsilon = \pi/6$.</figcaption></figure>

the spin component remains negligible. Figure 9 illustrates the variation of the CYORP obliquity component with the asteroid obliquity and the crater colatitude when $\delta$ and $\Delta$ 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.

# 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.
<figure><img src="image_9.png" /><figcaption>Fig. 9. Obliquity component (right panel) of the CYORP coefficient as a function of the asteroid obliquity $\epsilon$ and the crater colatitude $\alpha$ for $K = 2.65$ W m$^{-1}$ K$^{-1}$. Here, $d/D_0 = 0.16$ and $\delta = \Delta = 0$.</figcaption></figure>

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 Hanuš et al. (2013). To optimise the performance of our model, we assume a relatively high thermal conductivity of 1 W m$^{-1}$ K$^{-1}$. 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%.

## 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.} \quad (77) $$

On the other hand, the typical timescale associated with the YORP effect is approximately given by:

$$ \tau_{\text{YORP}} \sim 4 \left(\frac{R_{\text{ast}}}{1 \text{ km}}\right)^2 \left(\frac{\omega}{3 \times 10^{-4} \text{ s}^{-1}}\right) \text{Myr.} \quad (78) $$

Clearly, the YORP timescale is shorter than $\tau_{\text{imp,re-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_{\text{imp}}$ is

$$ \tau_{\text{imp}} = \frac{1}{P_i \pi R_{\text{Antonia}}^2 \Delta N(R > R_{\text{imp}})}, \quad (79) $$

where

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

Here, $P_i = 2.85 \times 10^{-18}$ km$^{-2}$ yr$^{-1}$ is the intrinsic collision probability, $C_R = 6 \times 10^5$, and $b_R = 2.2$ (Holsapple 2022). In the strength regime formulation, the crater produced by an impactor with the size of $R_{\text{imp}}$ has a size of

$$ R_{\text{crater}} = 1.306 R_{\text{imp}} \left(\frac{\rho_{\text{imp}}}{\rho_{\text{ast}}}\right)^{0.4} \left(\frac{Y}{\rho v_{\text{imp}}}\right)^{-0.205}, \quad (81) $$

with $Y = 100$ Pa and $v_{\text{imp}} = 5.3$ km s$^{-1}$. 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}, \quad (82) $$

where $g_{\text{ast}}$ is the surface gravity of the asteroid. In each time step (~$10^3$ 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
<figure><img src="image_10.png" /><figcaption>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.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>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.</figcaption></figure>

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}, \qquad (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{\epsilon} = \frac{T_{\epsilon}}{I\omega}, \qquad (84) $$

where $T_{\epsilon}$ 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 (Hanuš 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.

Union's Horizon 2020 research and innovation program under grant agreement no. 870377 (project NEO-MAPP).

7. Summary and conclusions

The YORP effect is a thermal torque produced by radiation
from the irregular surface of the asteroid. It has been demon-
strated that this effect is highly sensitive to surface topology
(Statler 2009; Breiter et al. 2009), including small-scale rough-
ness (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 conduc-
tivity. Using this model, we investigated the crater-induced
YORP (CYORP) effect in a computationally efficient man-
ner (about three orders of magnitude faster than the numeri-
cal method), allowing for a comprehensive exploration of the
functional dependence of the CYORP effect and its incorpora-
tion 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 W m⁻¹ K⁻¹).
This suggests that the model is suitable for application to materi-
als 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 inaccu-
racies. The typical value of the CYORP coefficient for the spin
component is 0.01, which is insensitive to the thermal param-
eter, 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 compu-
tation 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 coef-
ficient at each impact are assessed for the first time using our
CYORP model. While a complete investigation of the spin evo-
lution 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

References

Agarwal, J., Kim, Y., Jewitt, D., et al. 2020, A&A, 643, A152

Bottke, Jr, W. F., Vokrouhlický, D., Rubincam, D. P., & Nesvorný, D. 2006, Annu. Rev. Earth Planet. Sci., 34, 157

Bottke, Jr, W. F., Vokrouhlický, D., Walsh, K. J., et al. 2015, Icarus, 247, 191

Breiter, S., & Vokrouhlický, D. 2011, MNRAS, 410, 2807

Breiter, S., Bartczak, P., Czekaj, M., Oczujda, B., & Vokrouhlický, D. 2009, A&A, 507, 1073

Carruba, V., Domingos, R., Huaman, M., Santos, C. d., & Souami, D. 2014, MNRAS, 437, 2279

Carruba, V., Nesvorný, D., Aljbaae, S., & Huaman, M. E. 2015, MNRAS, 451, 244

Carruba, V., Nesvorný, D., & Vokrouhlický, D. 2016, AJ, 151, 164

Čuk, M., & Burns, J. A. 2005, Icarus, 176, 418

Čuk, M., Christou, A. A., & Hamilton, D. P. 2015, Icarus, 252, 339

Delbo, M., Walsh, K., Mueller, M., Harris, A. W., & Howell, E. S. 2011, Icarus, 212, 138

Đurech, J., Vokrouhlický, D., Pravec, P., et al. 2022, A&A, 657, A5

Farinella, P., Vokrouhlický, D., & Hartmann, W. K. 1998, Icarus, 132, 378

Fatka, P., Pravec, P., & Vokrouhlický, D. 2020, Icarus, 338, 113554

Golubov, O. 2017, AJ, 154, 238

Golubov, O., & Krugly, Y. N. 2012, ApJ, 752, L11

Golubov, O., & Lipatova, V. 2022, A&A, 666, A146

Golubov, O., & Scheeres, D. J. 2019, ApJ, 157, 105

Golubov, O., Scheeres, D., & Krugly, Y. N. 2014, ApJ, 794, 22

Golubov, O., Unukovych, V., & Scheeres, D. J. 2021, AJ, 162, 8

Hanuš, J., Đurech, J., Brož, M., et al. 2011, A&A, 530, A134

Hanuš, J., Đurech, J., Brož, M., et al. 2013, A&A, 551, A67

Holsapple, K. A. 2022, Planet. Space Sci., 219, 105529

Jacobson, S. A., & Scheeres, D. J. 2011, Icarus, 214, 161

Jacobson, S. A., Scheeres, D. J., & McMahon, J. 2013, ApJ, 780, 60

Jacobson, S. A., Marzari, F., Rossi, A., & Scheeres, D. J. 2016, Icarus, 277, 381

Lowry, S. C., Fitzsimmons, A., Pravec, P., et al. 2007, Science, 316, 272

Lowry, V. C., Vokrouhlický, D., Nesvorný, D., & Campins, H. 2020, ApJ, 160, 127

Lupishko, D., & Tielieusova, I. 2014, Meteor. Planet. Sci., 49, 80

Lupishko, D., Mikhalchenko, O., & Chiorny, V. 2019, Solar Syst. Res., 53, 208

Marzari, F., Rossi, A., Golubov, O., & Scheeres, D. J. 2020, AJ, 160, 128

Nesvorný, D., & Vokrouhlický, D. 2008, A&A, 480, 1

Polishook, D. 2014, Icarus, 241, 79

Pravec, P., Harris, A. W., Scheirich, P., et al. 2005, Icarus, 173, 108

Pravec, P., Harris, A., Vokrouhlický, D., et al. 2008, Icarus, 197, 497

Rozitis, B., & Green, S. F. 2012, MNRAS, 423, 367

Rozitis, B., & Green, S. F. 2013a, MNRAS, 433, 603

Rozitis, B., & Green, S. F. 2013b, MNRAS, 430, 1376

Rožek, A., Lowry, S., Nolan, M., et al. 2019, A&A, 631, A149

Rubincam, D. P. 2000, Icarus, 148, 2

Sánchez, P., & Scheeres, D. J. 2020, Icarus, 338, 113443

Scheeres, D. J. 2007, Icarus, 189, 370

Ševeček, P., Brož, M., Čapek, D., & Đurech, J. 2015, MNRAS, 450, 2104

Statler, T. S. 2009, Icarus, 202, 502

Steinberg, E. & Sari, R. 2011, AJ, 141, 55

Taylor, P. A., Margot, J.-L., Vokrouhlický, D., et al. 2007, Science, 316, 274

Tian, J., Zhao, H.-B., & Li, B. 2022, Res. Astron. Astrophys., 22, 125004

Veras, D., & Scheeres, D. J. 2020, MNRAS, 492, 2437

Vokrouhlický, D., & Čapek, D. 2002, Icarus, 159, 449

Vokrouhlický, D. & Nesvorný, D. 2008, AJ, 136, 280

Vokrouhlický, D., Milani, A., & Chesley, S. 2000, Icarus, 148, 118

Vokrouhlický, D., Nesvorný, D., & Bottke, W. F. 2003, Nature, 425, 147

Vokrouhlický, D., Brož, M., Bottke, W., Nesvorný, D., & Morbidelli, A. 2006, Icarus, 182, 118

Vokrouhlický, D., Bottke, W. F., Chesley, S. R., Scheeres, D. J., & Statler, T. S. 2015, Asteroids IV, eds. P. Michel, F. E. DeMeo, & W. F. Bottke (Tucson: University of Arizona Press), 509

Zegmott, T. J., Lowry, S., Rožek, A., et al. 2021, MNRAS, 507, 4914

Zhou, W.-H., Zhang, Y., Yan, X., & Michel, P. 2022, A&A, 668, A70