© The Authors 2026

1 Introduction

The scattering and re-emission of sunlight from the surface of an asymmetrically shaped asteroid may create a thermal torque, changing the asteroid's rotation period and spin axis direction in the long term. This phenomenon, known as the Yarkovsky-O'Keefe-Radzievskii-Paddack (YORP) effect (Rubincam 2000), has been adopted to explain some specific distributions of spin rates and axis orientations of asteroids (e.g. Pravec et al. 2008; Slivan et al. 2023). The YORP effect has already been directly confirmed by astronomical measurements on at least 12 asteroids (e.g. Taylor et al. 2007; Kaasalainen et al. 2007; Ďurech et al. 2008; Lowry et al. 2014; Hergenrother et al. 2019; Ďurech et al. 2022, 2024). However, contrary to the expectation of a 50:50 chance that an asteroid's spin rate may increase or decrease, all 12 samples exhibit positive angular velocity acceleration. Only recently, the asteroid (433) Eros is reported to rotationally decelerate (Feng et al. 2025). In fact, the numerical simulations also show that the YORP torque could increase or decrease the spin rate (e.g. Čapek & Vokrouhlický 2004), implying that the complicated process involved in the YORP effect may still not be fully understood.

Numerical simulations are often used to calculate the spin state evolution of asteroids driven by the YORP torque. Over the past decades, various numerical models have been developed based on different heat transfer models. Assuming zero thermal conductivity¹, Rubincam (2000) was the first to numerically model the YORP torque. Afterwards, using a 1D heat transfer model, in which heat conduction occurs only in the depth direction, Čapek & Vokrouhlický (2004) found the YORP-induced spin acceleration to be nearly independent of thermal inertia. This result was later confirmed by Breiter & Michalska (2008) and Nesvorný & Vokrouhlický (2008) using analytical models. Thus, when an asteroid's global shape is known, the 1D model can estimate the YORP effect fairly well and make accurate predictions (see e.g. Ďurech et al. 2008, for the case of (1862) Apollo).

Nevertheless, Statler (2009) and Breiter et al. (2009) suggested that small-scale topography could significantly influence the YORP torque through a shadowing effect. Rozitis & Green (2012) proposed the advanced thermophysical model, taking radiation blocking and re-absorption between surface elements into account. They found that the YORP torque depends sensitively on small-scale surface structures, even when their scale matches skin depth (typically a few centimetres for regolith).

In fact, at metre scales or smaller, the influence of asteroid surface topography becomes increasingly complex. Since numerical thermophysical models require a polyhedral shape model, the 1D assumption holds only when the skin depth is negligible compared to the surface facet size. Golubov & Krugly (2012) illustrated that decimetre-sized features (such as boulders), which have been neglected in previous studies, might produce a recoil force in the tangential direction, named the tangential YORP effect (TYORP). Since thermal waves can penetrate boulders, the 1D model is no longer valid at this scale. Ševeček et al. (2015) explained the anomalous acceleration of (25143) Itokawa by estimating the torque generated by boulders using a 3D heat transfer model. More recently, Nakano & Hirabayashi (2023) prove that lateral heat conduction could dampen or enhance the YORP torque depending on different shape models.

Impact craters are common small-scale features on asteroids, and their concave geometry naturally leads to complex radiative conditions. Recently, Zhou et al. (2022); Zhou & Michel (2024) find that even symmetric craters may generate a non-zero crater-induced YORP (CYORP) torque. This torque arises because the misalignment between the local normal and the radial direction of the crater, combined with shadowing and self-heating effects, breaks the symmetry of thermal re-emission from the crater surface, resulting in a net YORP torque.

In addition to the geometry-driven component investigated by Zhou et al. (2022); Zhou & Michel (2024), the YORP torque of a crater may also be attributed to internal heat conduction. The latter has received less attention but may become non-negligible when the skin depth is comparable to the crater size. In this study, we focused on metre-sized craters and investigated how lateral heat conduction contributes to the YORP effect induced by craters. To this end, we simulated the YORP effect created by craters using a complete thermophysical model, taking into account shadowing effect, self-heating effect, and lateral heat conduction.

The rest of this paper is organized as follows. In Sect. 2, we review the basic thermal model and introduce the numerical setting of our models. In Sect. 3, we investigate the dependence of YORP torque on the craters' locations, shapes, thermal properties, and rotation periods. The influence of lateral heat conduction is then explored in Sect. 4. Finally, we summarise our conclusions in Sect. 5.

2 Model and method

In this study, we used a toy model to investigate the YORP effect generated by a crater on a spherical asteroid. The YORP torque was calculated numerically after the temperature distribution on the surface was known, while the temperature was computed through a commercial numerical simulation software, COMSOL Multiphysics®².

2.1 YORP torque

In numerical models, the asteroid surface is generally treated as a polyhedron composed of triangular facets. For an illuminated facet i with area dS_i, photons re-emitting from it impart a recoil force df_i. Assuming a Lambert surface, the recoil force of the thermal radiation is given by $$\text{d}{f}_i=-\frac{2\in \sigma T_i^4}{3c}d{S}_i\cdot n_i,$$(1)

where ϵ is the thermal emissivity, σ is the Stefan-Boltzmann constant, T_i is the surface temperature, c is the speed of light, and n_i is the unit vector normal to the surface. The net torque T acting on the asteroid is the summation of all the facets: $$T={\displaystyle \sum _{}^N{r}_i\times }\text{\hspace{0.17em}d}{f}_i,$$(2)

where r_i is the position vector of facet i and N is the total number of facets. We adopted a widely accepted assumption that the body rotates about its principal axis of maximum moment of inertia. Thus, the torque vector can be decomposed into the spin component T_z and obliquity component T_ɛ. They introduce variations to the angular velocity ω and the obliquity ɛ respectively as $$\frac{\text{d}\omega }{\text{d}t}=\frac{T_z}{I_z}=\frac{T\cdot e}{I_z},$$(3)

and $$\frac{\text{d}\epsilon }{\text{d}t}=\frac{T_{\in }}{I_z\omega }=\frac{T\cdot e_1}{I_z\omega },$$(4)

where I_z is the moment of inertia around the spin axis (see e.g. Čapek & Vokrouhlický 2004). The vector e₁ is given by $$e_1=\frac{\left(N\cdot \text{e}\right)e-N}{\sin \epsilon },$$(5)

with N and e being the unit vectors normal to the orbit plane and the unit spin vector, respectively. Since the evolution timescale of the YORP effect is typically much longer than the asteroid's rotation and orbital periods, the torque in Eqs. (3) and (4) can be averaged over the rotation and orbital periods.

2.2 Basic model

For an asteroid spinning around its principal axis, the temperature T_i of the surface facets in Eq. (1) depends on the periodic external solar radiation input and the heat conduction process inside the body. We used COMSOL^® to simulate the heat transfer and thermal radiation in the fully 3D model. After obtaining the surface temperature distribution, the recoil force of thermal radiation was then calculated. The radiation interactions between the surface facets including shadowing, scattering, and re-absorption, as well as the lateral heat conduction beneath the surface were considered in the simulations. Distinct from the 1D model often used in previous studies, the 3D model can accurately describe the heat transfer process inside the body, particularly for the small-scale structures discussed in this paper.

We adopted a simplified model of a perfectly spherical asteroid of radius R = 10 m decorated with a single crater on its surface. Due to its perfect symmetry, a spherical surface produces no torque at all; thus, the total torque applied on the asteroid is simply the YORP torque generated by the crater. In addition, the alignment between the local normal of the spherical surface and the crater's radial vector ensures that the crater's orientation contributes no spin torque. This allows the effect of heat conduction to be studied independently of geometric contributions. The crater is assumed to be circular with a paraboloidal profile, for which the depth z as a function of radius r from its centre is given by $$z\left(r\right)=qD\left(1-\frac{4r^2}{D^2}\right)\times \{\begin{array}{ll}1,\hfill & \text{if\hspace{0.17em}}r\le D/2;\hfill \\ \exp \left[-\frac{{\left(r-D/2\right)}^2}{2{\delta }^2}\right],\hfill & \text{if\hspace{0.17em}}D/2<r\le D/2+3\delta ,\hfill \end{array}$$(6)

where D and δ are the diameter and the rim width, respectively (Statler 2009). The maximum depth is attained at the crater centre z_max = z(r = 0), and we define a dimensionless parameter, q = z_max/D, to represent the crater depth throughout this paper unless specified otherwise.

Under conditions of periodic incident radiation, heat conduction and temperature variations mainly occur at a very shallow depth beneath the surface, while the temperature remains constant deep inside (Wesselink 1948). Therefore, a crater-shaped shell with a certain thickness was adopted as the thermal model of the crater. The shell was set to be isothermal at the beginning of each simulation, with an initial temperature of 280 K determined after some test runs. The thermal conductivity κ = 0.001 W m⁻¹ K⁻¹, heat capacity C = 680 J kg⁻¹ K⁻¹, and surface density ρ = 1500 kg m⁻³ were assumed to be constant in our simulations. These are the typical values of regolith asteroids (Farinella et al. 1998). With a diameter of D = 2 m and a depth of q = 0.2, the crater is relatively small in size compared to the spherical asteroid with a 10 m radius. In fact, a simple algebraic calculation gives a ratio of the crater volume to the asteroid volume ~3q(D/R)³/32, which is 1.5 × 10⁻⁴ in this case.

The asteroid follows a circular orbit at a heliocentric distance of 1 AU, as it spins about its principal axis with a rotation period of P = 1000 s. We note that this rotation period is fairly short for asteroids with a 10 m radius.

In fact, for an asteroid with given thermal parameters (ρ, C, κ), its angular velocity ω = 2π/P determines a parameter Θ that measures relaxation between the absorption and re-radiation of energy (e.g. Lagerros 1996; Farinella et al. 1998; Xu et al. 2020): $$\Theta =\frac{\sqrt{\kappa \rho C\omega }}{\in \sigma T_{\text{s}s}^3},$$(7)

where T_ss is the sub-solar temperature. For a fast rotator with large ω, the relatively large Θ indicates a long relaxation time, resulting in small surface temperature variations and generally weak YORP torque. On the contrary, a high temperature gradient on the surface is generated when Θ is low. In this case, much higher mesh resolution is required to maintain accuracy in numerical simulations at significantly higher computational costs. Therefore, we adopted a short rotation period. As we show in Sect. 3, this period does not affect our qualitative results. This compromise between the model and the computational source has been applied in previous studies (e.g. Xu et al. 2022; Zhang et al. 2025).

2.3 Mesh model

The finite element method used in COMSOL requires a high-quality mesh model to ensure accurate simulation results. To find an appropriate mesh model, we ran test computations using a spherical asteroid model (without craters). As the asteroid rotates in sunlight, internal heat transfer occurs, and temperatures at different locations vary with time. Generally, the heat propagation distance can be characterised by the so-called 'skin depth' or 'penetration depth', $$l_s=\sqrt{\frac{\kappa }{\omega \rho C},}$$(8)

which represents the depth at which the amplitude of temperature variation decreases by a factor of e⁻¹ from the heat source (asteroid surface).

We assumed the asteroid has the same size and parameters as the asteroid introduced in Sect. 2.2, except that we also tested an additional thermal conductivity value of κ = 0.01 W m⁻¹ K⁻¹, alongside κ = 0.001 W m⁻¹ K⁻¹. Using a mesh of triangular-prisms with a layer thickness of l_s/2, we simulated the temperature evolution in the asteroid using COMSOL. We arbitrarily selected an equatorial point and show in Fig. the temperature variation along the radius through this point at sunset.

As shown in Fig. 1, the surface begins to cool down at sunset, with the temperature reaching a maximum at a certain depth below the surface. As temperature varies with depth, the temperature difference between adjacent mesh layers decreases rapidly with depth. The depths d₁ = 15l_s and d₂ = 12l_s for κ = 0.001 and 0.01 W m⁻¹ K⁻¹, respectively, where the temperature variation relative to the following deeper layer is less than 0.1%, are indicated by solid and dashed lines in Fig. 1. Below depths d₁ or d₂, the temperature remains nearly constant. Therefore, to simulate the heat transfer and temperature distribution in the asteroid, we used two types of meshes in our numerical model: a high-resolution mesh from the surface to a depth of 20l_s (> d₁, d₂) for accurate surface temperatures, and a low-resolution mesh in the interior to save on computational costs.

Specifically, the crater model adopted in this study consists of low-resolution tetrahedrons in the interior and a high-resolution 'shell' composed of triangular prisms beneath the surface, as sketched in Fig. 2. As we show in Fig. 1, the temperature is constant deep within the body, so we assume an isothermal core from a depth 100l_s. The bottom and side edges of the crater were set to be thermally adiabatic.

Furthermore, we consider heat transfer equilibrium reached when the temperature of any surface point varies by less than 0.1% between two consecutive passages through the subsolar point within one rotation period. In our simulation – which differs substantially from the 1D model – regions lacking illumination require more time to reach a steady temperature variation pattern due to lateral heat conduction. After some test runs, we find that equilibrium is reached in 5–20 rotation periods, depending on the thermal parameters.

Fig. 1 Temperature variations along a radius on the equatorial plane at sunset. Two thermal conductivities κ = 0.001 and 0.01 (given in W m−1 K−1) are adopted. The depth is measured downwards from the surface. The solid and dashed lines denote the depths from which the temperature variation relative to the subsequent mesh layer is less than 0.1% for κ = 0.001 and 0.01, respectively.

Fig. 2 Mesh model for the crater. The interior is composed of tetrahedrons, wrapped by a layer of triangular prisms near the surface.

Fig. 3 YORP accelerations (dω/dt, left panel; dɛ/dt, right panel) of the asteroid due to the crater as a function of obliquity ɛ. The different colours show the crater at ten different latitudes from ϕ = 0° to ϕ = 90°, as indicated by the colour bar.

3 Spin and obliquity torques

In the simple model of a spherical asteroid decorated with a crater, the YORP effect arises solely from the torque produced by the crater. In addition, since the crater's volume and mass are negligible compared to those of the asteroid (the ratio between them is ~ 10⁻⁴), we assume that the changes in the moment of inertia and principal axis due to the crater are both negligible.

Using COMSOL, we simulated the thermal evolution of the asteroid with a crater from an initial isothermal state. After tens of rotation periods, when a dynamic equilibrium was attained inside the body, we calculated the recoil force following Eq. (1). Knowing the recoil force, the torque may be computed and then averaged over the rotation and revolution periods. The asteroid's angular velocity ω and obliquity e may be altered by this torque, and the rates of change (dω/dt, dɛ/dt) may be used to measure the strength of the YORP effect.

3.1 Crater at different locations

For an asteroid spinning about an axis tilted by ɛ relative to the orbital plane normal, the illumination conditions of a crater vary with its location on the surface, as does the YORP torque it generates. Since the net effective torque is averaged over the rotation period, only the latitude of the crater matters.

Setting a crater with D = 2 m, q = 0.2 at different latitudes on the surface of an asteroid with R = 10 m, we simulated the thermal dynamics and obtained the YORP torque, from which dω/dt and dɛ/dt were calculated. We summarise our results in Fig. 3.

As shown in the left panel of Fig. 3, dω/dt ≥ 0 at all latitudes ϕ, indicating that the crater always accelerates the asteroid's spin rate wherever it locates. From ϕ = 90° with the crater at the pole where dω/dt = 0, the YORP torque increases monotonically to ϕ = 0° when the crater lies on the equator. Since the spin torque component depends on incident energy (e.g. Čapek & Vokrouhlický 2004; Nesvorný & Vokrouhlický 2008; Golubov et al. 2016), the spin acceleration weakens when at higher latitudes where the crater receives less solar radiation due to lower incident angles and stronger self-shadowing. We note that the maximum dω/dt at ϕ = 0° also implies that the crater in the south hemisphere will accelerate the spin rate.

For a crater at a given latitude, the spin acceleration decreases as obliquity ɛ increases, reaching a minimum at ɛ = 90° when the asteroid Ties' on the orbital plane. Beyond this, the variation reverses for 90° ≤ ɛ ≤ 180°. We note some unexpected variations in our calculations, particularly for ϕ = 20° near ɛ = 0° (180°), which we attribute to calculation errors.

The right panel of Fig. 3 shows the crater-induced variation in obliquity. In contrast to the reflectional symmetry about ɛ = 90° seen in the left panel, the curves in the right panel exhibit rotational symmetry with respect to ɛ = 90°. For a pro-grade spinning asteroid with a given obliquity (ɛ < 90°), dɛ/dt changes from negative to positive as the crater moves from the equator to the north pole, indicating that polar craters generally increase obliquity while equatorial craters decrease it.

For the asteroid with a crater at a specific location (here indicated by its latitude ϕ), the YORP effect may drive its obliquity to an equilibrium state, where dɛ/dt = 0. As shown in the right panel of Fig. 3, when the crater is near the equator (ϕ = 0°, 10°, 20°), ɛ = 0° is a stable equilibrium or 'asymptotic obliquity' (Vokrouhlický & Čapek 2002), since dɛ/dt < 0 for ɛ > 0° around this point. For craters at high latitudes (ϕ ≥ 50°), however, ɛ = 0° becomes an unstable equilibrium, since dɛ/dt > 0 for ɛ > 0°.

Conversely, ɛ = 90° is an unstable equilibrium for near-equatorial craters but an asymptotic obliquity for high-latitude craters. For mid-latitude craters (in Fig. 3, ϕ = 30°, 40°), the dɛ/dt curve crosses zero at ɛ ~ 67° and ~38° for ϕ = 30° and 40°, respectively. This additional equilibrium is unstable, and its appearance makes both ɛ = 0° and ɛ = 90° asymptotic obliquities. Namely, when the crater is at mid-latitudes (e.g. ϕ = 30°), the YORP torque tilts the spin axis towards ɛ = 90° if the initial obliquity ɛ₀ ≳ 67° or towards ɛ = 0° if ɛ₀ ≲ 67°.

In summary, a low-latitude crater produces a YORP torque that drives the spin axis perpendicular to the orbital plane (ɛ = 0°), while a high-latitude crater tilts the spin axis towards the orbital plane (ɛ = 90°). For a mid-latitude (~30°–50°) crater, it yields two asymptotic obliquities, either 0° or 90°, depending on its initial state.

Note that the spin torque from a single crater is weak. Specifically, for ɛ = 0° and ϕ = 0°, the equatorial crater produces a spin torque of $T_z^{\epsilon ,\varphi =0}~1.2\times {10}^{-8}\text{N\hspace{0.17em}m}$ Nm. Crater torques at any latitude can be estimated by interpolating the nine values used in Fig. 3. Assuming 200 such craters distributed randomly across the asteroid, covering half of the total surface area, their combined spin torque would be about 100 times that of $T_z^{\epsilon ,\varphi =0}$, i.e. $T_z^{\epsilon ,\varphi =0}~1.2\times {10}^{-6}\text{N\hspace{0.17em}m}$. Golubov et al. (2014) suggested normalizing the YORP torque T_z via $${\tau }_z=\frac{T_{z}c}{\Phi r_{\text{eq}}^3},$$(9)

where c, Φ, and r_eq denote the light speed, the solar radiation flux, and the equivalent radius of asteroid, respectively. To distinguish boulder-induced TYORP, Golubov et al. (2014) named the torque from the same terrain without boulders 'normal YORP' (NYORP). The dimensionless normalized spin torque of our half-crater-covered asteroid is τ_z ~ 2.6 × 10⁻⁴. This value falls within the typical NYORP range for real asteroids, which is of the order of ~10⁻⁴ to 10⁻³ (Golubov et al. 2014). Thus, for small asteroids such as those in our model, the crater-induced spin torque is comparable to NYORP, roughly 10–100% of typical NYORP values. It is important to note that this spin torque is consistently positive, definitely accelerating the asteroid's rotation.

Conversely, such randomly distributed craters are unlikely to produce net obliquity YORP torque because the torques from the northern and southern hemispheres typically cancel out, resulting in little to no net effect.

Fig. 4 As in Fig. 3, but for different crater depths q (as indicated by the colour bar) and crater location at latitude ϕ = 40°.

3.2 Crater of different shapes

For a crater of a given diameter D and profile defined by Eq. (6), its shape is characterized by the depth q. As q varies, the crater surface morphology changes, and consequently, the incident and emitted radiations vary accordingly. Below, we discuss how the YORP torque varies with crater depth.

The same model was adopted: a crater of diameter D = 2 m on an asteroid of radius R = 10 m. The location of the crater was fixed at ϕ = 40°. A crater disappears when its depth approaches zero (q → 0), while craters with q > 0.5 (a hemispherical crater has q = 0.5) are very unlikely to exist. Therefore, besides q = 0.2 discussed in Sect. 3.1, we repeated the calculations of dω/dt and dɛ/dt for four other depth values q = 0.1, 0.3, 0.4, and 0.5 and summarise the results in Fig. 4.

Figure 4 illustrates that the spin acceleration dω/dt and obliquity change rate dɛ/dt both positively correlate with the depth q. Specifically, increasing q from 0.1 to 0.5 significantly enhances both by approximately one order of magnitude. A higher q value indicates a deeper crater with a relatively steeper wall. And the same amount of irradiation energy absorbed and then emitted by a steep crater produces a stronger torque than a shallow crater.

As shown in Fig. 3, a mid-latitude crater gives rise to an additional unstable equilibrium of obliquity. In the right panel of Fig. 4, a larger crater depth drives the unstable equilibrium towards ɛ = 0° (or 180°). As a result, the asymptotic obliquity of 90° is more likely and more quickly to occur for asteroids with deep craters.

For a crater at a specific location, the diameter determines the solar-illuminated and thermal-radiating areas. Therefore, the YORP torque component arising from the overall thermal radiation along the crater's orientation – that is, the obliquity component (dɛ/dt) – is simply proportional to the square of the diameter. But for the YORP torque (dω/dt) arising from symmetry breaking due to internal heat conduction (see Sect. 4), a smaller diameter implies greater asymmetry within the crater, although the total radiation is correspondingly smaller. This makes the influence of the crater size on the YORP effect complicated.

3.3 Thermal conductivity

In a 1D heat transfer model, heat conducts only along the depth direction, penetrating only a few skin depths, l_s, which is generally much smaller than the size of the asteroid. For each surface element on the crater in a 1D model, it absorbs solar radiation, transfers energy downwards, and then after a certain delay emits the same amount of energy through irradiation from the same surface element. Considering crater symmetry, the net spin torque (T_z) produced by the crater in such a 1D thermal model is in fact null. But in a 3D model, the heat diffuses in all possible directions, and the spin torque leads to an acceleration dω/dt, as our numerical results (Figs. 3 and 4) have shown.

The thermal inertia $\left(\Gamma =\sqrt{\kappa \rho C}\right)$ characterises an asteroid's thermal response, where a high Γ reduces temperature variations and introduces a significant 'thermal lag' by delaying the thermal response to radiation. As the resulting heat transport dictates the net YORP torque, the torque produced by a crater is expected to depend on Γ. To quantify this, we performed additional calculations for a crater at 40° latitude, varying Γ by changing κ between 0.001 and 0.01 W m⁻¹ K⁻¹ while keeping ρ and C fixed. The results in Fig. 5 illustrate this dependence.

Higher thermal conductivity facilitates faster heat diffusion in all directions, thereby creating greater deviation from the 1D model where the net spin torque is zero. Consequently, as shown in the left panel of Fig. 5, the spin acceleration dω/dt increases with κ. Notably, the differences between the adjacent curves dimmish as κ rises, suggesting that the spin acceleration approaches an upper limit. This saturation is consistent with the theoretical expectation for infinite thermal inertia (Γ → ∞), where temperature variations vanish and the YORP torque drops to zero. Our computations, however, were not extended to this high Γ regime due to limited computing resources.

As shown in the right panel of Fig. 5, the influence of κ on the obliquity component of the YORP effect is similar to the influence of crater depth seen in Fig. 4: higher κ values lead to higher probabilities of asymptotic obliquity at 90°. However, it is worth noting that the obliquity YORP torque exhibits considerably less sensitivity to thermal conductivity. This is evidenced by the markedly smaller variation in dɛ/dt with thermal conductivity in Fig. 5 compared to Fig. 4. This insensitivity is further supported by the much smaller relative change in dɛ/dt than in dω/dt in Fig. 5.

Fig. 5 As in Fig. 3, but for different thermal conductivities (as indicated by the colour bar). The crater has q = 0.2 and is located at ϕ = 40°.

Fig. 6 Spin YORP acceleration produced by an equatorial crater as a function of rotation period P. The other model parameters are the same as in Fig. 3.

3.4 Rotation period

Up to this point, a short asteroid rotation period P = 1000 s was adopted, as a compromise between the model and computation cost. As a matter of fact, the steady-state rotation periods of rubble-pile asteroids typically fall within the range of hours. Keeping other parameters of the model as before but varying the rotation period in 2000–20 000 s, we calculated the spin YORP torque of an equatorial crater. As summarised in Fig. 6, the spin torque decreases with rotation period. For P = 20000 s, the dω/dt is nearly one order of magnitude smaller than that for P = 1000 s.

According to the thermophysical model (Lagerros 1996), surface temperature is determined by the thermal parameter 0 once the shape is fixed, while 0 represents the ability to maintain the temperature during periodic external insolation (Harris & Drube 2020). When the rotation period is fairly long for large asteroids, their Θ ≈ 0, and the spin YORP effect is dominated by geometric asymmetry rather than internal heat conduction, as in the case here.

4 Lateral heat conduction

For any structure on an asteroid's surface, because the lever arm and normal direction of any surface element dS remain constant in a body-fixed frame, the rotation-averaged torque produced by dS solely depends on the re-emitted energy, which in the 1D model is nearly equal to the total radiation received by the same surface element. In this scenario, the spin YORP torque is nearly independent of the thermal inertia (Nesvorný & Vokrouhlický 2008), and the symmetric topology of the crater in our model does not generate the spin YORP torque.

However, lateral heat conduction in the 3D model might break the balance of the energy budget. Consequently, our numerical results demonstrate positive spin acceleration dω/dt (Fig. 3) and its dependence on the thermal inertia (Fig. 5). Lateral heat conduction inside the body plays a critical role in creating the spin torque, as we show below.

4.1 Asymmetry in temperature distribution and YORP torque

For any surface element on the crater, the perfect equality between the absorbed solar radiation and the re-emitted thermal radiation could break if either lateral heat conduction or self-heating occurs. Particularly, a spin YORP torque may arise from the asymmetries in both the spatial temperature distribution and its temporal variation between the eastern and western sides of the crater. To show this, we located a crater of diameter D = 2 m and depth parameter q = 0.2 on the equator of an asteroid with radius 10 m. The north and south halves of such a crater are absolutely symmetrical to each other, so there is no obliquity YORP torque.

Arbitrarily, we selected two point pairs on the crater's prime vertical circle (the great circle passing through the east, nadir, and west points): one at (r = D/2,z = 0) and the other at (r = D/4,z = 30/20). According to Eq. (6), the first pair lies on the crater's edge where it meets the rim, while the second is located lower on the crater wall at a depth of z = 30 cm. The temperatures at these four points were calculated in the 3D model and plotted over one rotational period in Fig. 7, with the rotational phase starting from midnight.

During the local night, temperatures at all four points decline slowly under radiative cooling, then surge abruptly at morning illumination. The west-side points, receiving sunlight first at a high incidence angle, warm rapidly, triggering lateral subsurface heat conduction towards the east. As a result, the east-side points experience pre-sunrise warming due to both this heat transfer and the self-heating effect. At local sunrise, the east-side temperatures rise more gradually owing to the lower solar incidence angle, taking longer to peak. Additionally, the shadowing effect causes the east-side points to lose direct solar illumination at a higher incidence angle, leading to an abrupt temperature drop around the rotational phase of 3π/2. As summarised in Fig. 7, this asymmetry creates not only a time lag but also a pronounced difference in the temperature profiles between the two sides.

For symmetrically placed points on the east and west sides, the west point exceeds the east in temperature for only about a quarter of the rotational period. Furthermore, the east point exhibits higher minimum and maximum temperatures than its western counterpart. For instance, as shown in Fig. 7 for points at the rim (z = 0), the west point's temperature ranges from 203.6 to 381.0 K, whereas the east point's varies from 207.5 to 384.5 K.

A higher surface temperature leads to a stronger thermal radiation recoil force. This results in an imbalance where the positive spin YORP torque generated by the eastern half may exceed the negative torque from the western half. To quantify this, we computed the spin YORP torques from the two sectors separated by the crater's meridian. Their variations over a rotational period are shown in Fig. 8. Since this torque asymmetry originates partly from lateral heat conduction, an effect enhanced under higher thermal conductivity κ, we present calculations for two κ values for comparison.

Figure 8 clearly reveals a time lag between the positive spin torque (from the east half) and the negative torque (from the west half). A more subtle but crucial feature is the asymmetry in their temporal profiles. Compared to the case of κ = 0.001 W m⁻¹ K⁻¹, a higher thermal conductivity κ = 0.01 W m⁻¹ K⁻¹ (dashed lines) strengthens internal heat conduction, dampening the surface temperature variation. This reduces the peak torque magnitudes for both crater halves. Furthermore, enhanced conductivity allows more heat to be stored internally and released as thermal radiation during the night (approximately during phases (0, π/2) and (3π/2, 2π)). This effect produces a distinctive 'tail' in the net positive spin torque (red lines) after sunset (phase 3π/2), which constitutes the primary contribution to the total averaged YORP torque. Of course, it is worth noting that the torque might disappear also in the extreme (but highly improbable) case of extremely high thermal conductivity or a very small asteroid, where temperature variations would be negligible.

Fig. 7 Temperature variations at symmetrical points on the crater's prime vertical circle. The solid and dashed lines correspond to points at depths of z = 0 and z = 30 cm, respectively, while the blue and yellow colours indicate locations on the west and east sides of the crater.

Fig. 8 Spin YORP torques produced by the western (yellow) and eastern (blue) halves of the crater. The total torque (the summation of both halves) is also plotted in red. The solid and dashed lines represent thermal conductivity κ = 0.001 and 0.01 W m−1 K−1, respectively.

4.2 A pseudo-crater model

To elucidate how a positive YORP torque can arise from a symmetric structure such as a crater, we introduce a 'pseudo-crater' model. This model replaces the crater in Fig. 2 with a semi-cylindrical shell of radius r = 0.1 m, oriented north-south on the equator of a spherical asteroid with radius R = 10 m. As the asteroid rotates, every point along a cylindrical generatrix of this shell experiences identical solar radiation conditions. Consequently, no temperature gradient is established in this longitudinal direction, and thus, no heat diffusion occurs along it. Effectively, this reduces the system to a 2D model where the heat diffusion is confined to the east-west direction. For simplicity, we assigned the shell a unit length (h= 1 m) and the same thermal parameters as the crater model used previously in our study. Thermal diffusion and radiation from the cylinder's top and bottom ends were neglected.

To evaluate the YORP torque generated by the pseudo-crater, we computed the surface temperature distribution starting from an isothermal state. Subsequently, the spin YORP torque was obtained by averaging the recoil force torque over an entire asteroid rotation. We employed three distinct thermal models: a 1D model confined to vertical heat conduction; a 2D model incorporating lateral conduction (only in the east-west direction) but neglecting self-heating; and a full 2D model accounting for both lateral conduction and self-heating. A comparison of the resulting torques (Fig. 9) may help us isolate the contribution of lateral heat diffusion and self-heating.

Starting from an isothermal state, both the eastern and western facets lose energy through thermal radiation before sunrise. However, the eastern facet – which receives sunlight later – cools more and thus begins the first rotation at a lower temperature. Consequently, it stores more solar energy while emitting less thermal radiation than the western facet, resulting in an initial negative spin torque in both the 1D and 2D models. From the second rotation onwards, self-heating pre-warms the eastern facet before sunrise. This process transfers extra energy that must be radiated away, generating a positive net spin torque. In contrast, models without self-heating lack this energy transfer; as Fig. 9 shows, the net torque gradually converges to a certain value as heat diffuses inside the body towards a dynamic equilibrium.

After the initial rotation, the self-heating effects on the eastern and western crater facets become largely symmetric as the temperature difference between two sides diminishes. The net torque is thereby dominated by direct solar irradiation and internal thermal conduction. This is confirmed by the convergence of the time-averaged torques to stable values at dynamic equilibrium (Fig. 9). The close agreement between the final values in the 2D model without self-heating (4.75 × 10⁻¹⁰ N m) and the full 2D model (5.12 × 10⁻¹⁰ N m) strongly suggests that self-heating is inconsequential for the overall YORP torque in this scenario.

In contrast to the 2D models, the 1D model yields a final YORP torque of only 5.84 × 10⁻¹¹ N m, nearly an order of magnitude smaller and effectively negligible. A separate 1D simulation without self-heating (not shown) converges to a nearly identical value, confirming that the pseudo-crater generates only negligible spin torque when heat diffusion is restricted to the depth direction. The dramatic increase in torque observed in the 2D models unequivocally indicates lateral heat transfer as the primary driver of the spin YORP effect in this crater structure.

Since numerical experiments demonstrate that the spin YORP torque arises mainly from lateral heat transfer within the body, this torque varies with thermal conductivity κ in a given model. Further tests with different κ values in the same pseudo crater model revealed a notable dependence of the spin YORP torque on this parameter. In fact, when κ is very small, the asymmetric temperature distribution caused by thermal lag is negligible. On the other hand, the temperature difference between the west and east sides vanishes when κ is very large. Therefore, the maximal spin YORP torque is attained at some specific κ in between, which is found to be κ ~ 2 × 10⁻² W m⁻¹ K⁻¹, corresponding to a skin depth of l_s ≈ 0.2 cm, in our preliminary calculations. We leave a detailed investigation on the dependence of this torque on thermal parameters to a separate paper.

Fig. 9 Convergence of the calculated YORP torque produced by the pseudo-crater. The iteration number indicates the number of asteroid rotation periods elapsed since the beginning of our calculation (assuming an isothermal state).

5 Conclusions

While the YORP effect is known to depend on an asteroid's macroscopic shape and thermal parameters, its sensitivity to small-scale surface features, including boulders, craters, and roughness, is also recognized (e.g. Breiter et al. 2009; Statler 2009; Golubov & Krugly 2012; Rozitis & Green 2012). Recent studies of the CYORP effect showed that the YORP torque may stem from the crater's orientational asymmetry (Zhou et al. 2022; Zhou & Michel 2024). In these studies, lateral heat conduction inside the asteroid is often ignored. Thus, a fully symmetrical configuration – for example a rotationally symmetrical crater located on the equator of a spherical asteroid – produces no net torque. However, when internal heat conduction (an ever-present factor in reality) is considered, a net YORP torque can emerge. In this study, we employed a complete 3D model to calculate the radiation interactions between surfaces and heat conduction within the body. We demonstrated that heat conduction can indeed generate a net YORP torque, even in perfectly symmetrical crater geometries.

To eliminate the influence of geometric characteristics, we adopted a simple spherical asteroid model featuring a circular crater with a paraboloidal profile bounded by a raised rim. The radii of the asteroid and crater were set to R = 10 m and r = 1 m, respectively, while the crater morphology was characterised by its depth-to-diameter ratio q. Under the assumption of zero thermal inertia, this symmetrical configuration generates a torque in obliquity but yields no net spin torque over a full rotation.

Using typical thermal parameters, we numerically solved the 3D heat conduction equation to obtain the temperature distribution on the crater's surface. From this, the recoil force and YORP torque were calculated. Our results show that the spin torque generated by the crater is consistently positive, with its magnitude increasing monotonically as the crater's location shifts from the polar to the equatorial region (Fig. 3). The obliquity YORP torque, however, demonstrates a non-monotonic dependence on crater location. A low-latitude crater drives the obliquity towards an asymptotic equilibrium at ɛ = 0°. In contrast, a high-latitude crater drives it towards ɛ = 90°. At mid-latitudes (for example, ϕ = 30° and 40°), both 0° and 90° are stable equilibrium points for obliquity e (Fig. 3).

The crater's geometry also affects the torque. As the depth parameter q increases from 0.1 to 0.5, changing illumination conditions within the crater alter both the spin and obliquity torques by one order of magnitude. Specifically, a larger q indicates steeper crater walls, leading to a larger temperature gradient on the surface. Both heat conduction within the body and the shadowing effect become more pronounced at higher q, and consequently, both the spin and obliquity torques increase (Fig. 4).

In our model of a spherical asteroid decorated with a symmetrical crater, the obliquity torque is primarily produced by the overall thermal radiation from the crater, determined mainly by the crater's orientation (aligned with the position direction) and the shadowing effect. The strength of these obliquity torques, measured in normalized torque, ranges from ~10⁻⁶ to 10⁻³, consistent with the findings of Zhou et al. (2022). However, the spin torque, which is about two orders of magnitude smaller than the obliquity torque (Fig. 3), arises mainly from symmetry breaking due to lateral heat conduction beneath the crater's surface. Therefore, an increase in thermal conductivity significantly enhances the spin torque, while the obliquity torque is only marginally affected (Fig. 5).

For the same reason, crater size (diameter) directly affects the obliquity torque: it is approximately proportional to the square of the diameter. The dependence of the spin torque on diameter, however, is somewhat more complex. For sufficiently large craters, the asymmetry in surface temperature distribution caused by lateral heat conduction approaches a limiting state such that the spin torque tends to an asymptotic value $M_z^{\infty }$. Our preliminary calculations indicate that the spin torque M_z depends on the crater diameter D, as $M_z\approx M_z^{\infty }/\left[1+{\left(D_0/D\right)}^{3/2}\right]$. A detailed investigation is left for future work.

In our model, we set a rotation period P = 1000 s. For longer periods, the results qualitatively remain the same. Quantitatively, the spin YORP torque decreases to about one-sixth as P increases from 1000 s to 20 000 s (Fig. 6).

A crater on a spherical asteroid equator, with its orientation aligned with the local surface normal, produces no obliquity torque; nevertheless, it serves as a clear illustration of spin torque generation. We calculated temperatures at representative points in such a crater (Fig. 7) and the resulting spin torques from thermal radiation on its eastern and western halves (Fig. 8). The temporal profiles of these temperatures and torques reveal how an asymmetric temperature distribution develops in this symmetrical crater configuration, ultimately producing a net positive spin torque.

To isolate the essential influence of lateral heat diffusion on spin YORP torque generation, we constructed a pseudo-crater model where the temperature gradient was confined to the east-west direction. Our numerical simulations demonstrate that spin torque arises exclusively when lateral heat conduction is allowed (in a 2D model), with the self-heating effect making only a negligible contribution (Fig. 9).

We quantified the spin YORP torque arising from lateral heat conduction in a crater using the dimensionless normalized torque and find it to be comparable to the normal YORP. Consequently, the cumulative effect from multiple craters is non-negligible. A rough estimate suggests that for a heavily cratered asteroid, the cumulative contribution from craters could range from 10 to 100% of the normal YORP torque.

In our calculations, crater orientation effect was not considered. We assumed crater orientation to be aligned with the local surface normal. It is worth noting that the obliquity component of the YORP torque produced by craters located off the equator (right panels of Figs. 3, 4 and 5) is, to some extent, physically equivalent to the torque generated by craters with different orientations, as studied in Zhou et al. (2022).

We note that the calculations in this paper are mainly based on an ideal model, in which a crater of radius 1 m sits on an asteroid of radius 10 m. The lateral-heat-conduction-induced spin YORP torque might be suppressed by geometric effects for large craters and asteroids. The YORP effect on real asteroid is often more complex due to the combined contribution of TYORP, CYORP, NYORP, and the boulder-induced YORP (Baker & McMahon 2025). The general solution of geometric and thermophysical effects across scales therefore stands as a key priority for future work in this field.

Acknowledgements

This work has been supported by the National Natural Science Foundation of China (NSFC, Grants No.12373081 & No.12150009) and the China Manned Space Program with grant No.CMS-CSST-2025-A16.

References

All Figures

	Fig. 1 Temperature variations along a radius on the equatorial plane at sunset. Two thermal conductivities κ = 0.001 and 0.01 (given in W m−1 K−1) are adopted. The depth is measured downwards from the surface. The solid and dashed lines denote the depths from which the temperature variation relative to the subsequent mesh layer is less than 0.1% for κ = 0.001 and 0.01, respectively.
In the text

	Fig. 2 Mesh model for the crater. The interior is composed of tetrahedrons, wrapped by a layer of triangular prisms near the surface.
In the text

	Fig. 3 YORP accelerations (dω/dt, left panel; dɛ/dt, right panel) of the asteroid due to the crater as a function of obliquity ɛ. The different colours show the crater at ten different latitudes from ϕ = 0° to ϕ = 90°, as indicated by the colour bar.
In the text

	Fig. 4 As in Fig. 3, but for different crater depths q (as indicated by the colour bar) and crater location at latitude ϕ = 40°.
In the text

	Fig. 5 As in Fig. 3, but for different thermal conductivities (as indicated by the colour bar). The crater has q = 0.2 and is located at ϕ = 40°.
In the text

	Fig. 6 Spin YORP acceleration produced by an equatorial crater as a function of rotation period P. The other model parameters are the same as in Fig. 3.
In the text

	Fig. 7 Temperature variations at symmetrical points on the crater's prime vertical circle. The solid and dashed lines correspond to points at depths of z = 0 and z = 30 cm, respectively, while the blue and yellow colours indicate locations on the west and east sides of the crater.
In the text

	Fig. 8 Spin YORP torques produced by the western (yellow) and eastern (blue) halves of the crater. The total torque (the summation of both halves) is also plotted in red. The solid and dashed lines represent thermal conductivity κ = 0.001 and 0.01 W m−1 K−1, respectively.
In the text

	Fig. 9 Convergence of the calculated YORP torque produced by the pseudo-crater. The iteration number indicates the number of asteroid rotation periods elapsed since the beginning of our calculation (assuming an isothermal state).
In the text