© ESO, 2013

1. Introduction

The effect of sunlight on the dynamics of small asteroids is now known to be considerable over long time spans. The anisotropic reflection and thermal re-emission of solar photons is essentially equivalent to a large number of small rocket engines distributed on the surface of an asteroid. Thus the dynamical mechanisms are usually divided into two categories: the net propulsion force changing the orbital motion (the Yarkovsky effect), and the net torque affecting the rotational dynamics (the Yarkovsky-O’Keefe-Radzievskii-Paddack or YORP effect). Both have been observed directly (Chesley 2003; Kaasalainen et al. 2007; Lowry et al. 2007; Ďurech et al. 2008, 2012), and there is clear indirect evidence of their long-term role in the evolution of asteroid orbits and spin states (Bottke et al. 2001; Vokrouhlický et al. 2003, 2006).

The analytical modeling of the YORP torques, or any other illumination-dependent functions, on general shapes is constrained by the concept of insolation: in which geometries is a surface patch visible? For convex bodies insolation is trivial, and they are the only ones for which exact analytical torques are possible. Various semianalytical approaches have developed small nonconvex perturbations while neglecting numerical ray-tracing within some perturbation size (e.g., Breiter & Michalska 2008; Nesvorný & Vokrouhlický 2007, 2008a; hereafter NV07 and NV08). Unfortunately, the shadowing error cannot be estimated in terms of the perturbation parameter, so no validity regime can be established for such approximations even when the perturbation is very small.

We present the fundamentals of the YORP effect and its stability and symmetry properties in a conceptually and technically simple analytical and exact formulation. For example, it is easy to identify the enigmatic YORP instability with the “trash coefficients” of a Laplace series, and to show that, while quadrant-symmetric bodies have no secular YORP torques for zero thermal inertia (K = 0), they acquire a nonvanishing component for K ≠ 0. Breiter et al. (2011; hereafter BRV11) arrive at some of the results of Sect. 3 (K = 0) by a different context and approach. We provide easy-to-use formulae directly applicable to typical body representations such as polyhedra. The motion-averaged formulation pertains to any shapes.

We discuss the basic concepts and functions in Sect. 2. The cases of zero and nonzero thermal inertia are studied in, respectively, Sects. 3 and 4. The YORP effect can be very sensitive to the details of the shape and other physical characteristics of the body (Statler 2009; Rozitis & Green 2012), so we discuss YORP stability in Sect. 5. We sum up in Sect. 6, and numerical formulae for convex polyhedral shape representations are presented in an appendix.

2. Problem setup and definitions

Let the direction of the illumination source in a coordinate system fixed to the target body (the z-axis coinciding with the rotation axis) be denoted by ω ∈ S². Here entities on the unit sphere S², defined by two direction angles, are identified with unit vectors in R³. Thus, e.g., the outward unit normal vector η ∈ S² is given by η = η(ϑ,ψ) (with ϑ measured from the pole), 0 ≤ ϑ ≤ π, 0 ≤ ϑ < 2π; i.e., ${\mathit{\eta }}_{\mathrm{1}}\mathrm{=}\sin \mathit{\vartheta }\cos \mathit{\psi ,} {\mathit{\eta }}_{\mathrm{2}}\mathrm{=}\sin \mathit{\vartheta }\sin \mathit{\psi ,} {\mathit{\eta }}_{\mathrm{3}}\mathrm{=}\cos \mathit{\vartheta }\mathit{.}$(1)Likewise, ω = ω(θ,ϕ) (the subsolar coordinates). The position of the asteroid around the Sun in the orbital frame is given by λ, 0 ≤ λ < 2π, and the tilt (obliquity) of the rotation axis from the orbital pole is denoted by ϵ, 0 ≤ ϵ ≤ π. The orbital frame and the pole direction of the asteroid are chosen such that $\cos \mathit{\theta }\mathrm{=}\mathrm{\left(}\cos \mathit{\lambda ,}\sin \mathit{\lambda ,}\mathrm{0}\mathrm{\right)}\mathrm{·}\mathrm{\left(}\sin \mathit{ϵ,}\mathrm{0}\mathit{,}\cos \mathit{ϵ}\mathrm{\right)}\mathrm{=}\cos \mathit{\lambda }\sin \mathit{ϵ}\mathit{.}$(2)For clarity, we assume a circular orbit around the Sun. Ellipticity is straightforward to include as described in NV07 and NV08, and in any case the modifications due to ellipticity are only of order $\mathrm{(}\sqrt{\mathrm{1}\mathrm{-}{\mathit{e}}^{\mathrm{2}}}{\mathrm{)}}^{-1}$ for eccentricity e (NV07). Denoting the rotation angle of the asteroid around its axis by ϕ′: = Ωt, the transformation of a vector v from the orbital frame to the body frame reads ${{v}}_{\mathrm{ast}}\mathrm{=}{{R}}_{\mathit{z}}\mathrm{\left(}{\mathit{\varphi }}^{\mathrm{\prime }}\mathrm{\right)}{{R}}_{\mathit{y}}\mathrm{\left(}\mathit{ϵ}\mathrm{\right)}{{v}}_{\mathrm{orb}}\mathit{,}$(3)where R_i are the usal rotation matrices for z- and y-axes (Kaasalainen & Lamberg 2006). The transformations to an oblique frame aligned with the rotation axis but not rotating with the asteroid are ${{v}}_{\mathit{ϵ}}\mathrm{=}{{R}}_{\mathit{z}}\mathrm{(}\mathrm{-}{\mathit{\varphi }}^{\mathrm{\prime }}\mathrm{)}{{v}}_{\mathrm{ast}}\mathrm{=}{{R}}_{\mathit{y}}\mathrm{\left(}\mathit{ϵ}\mathrm{\right)}{{v}}_{\mathrm{orb}}\mathit{.}$(4)The oblique frame is needed in computing the torque component that changes the obliquity of the spin axis. We assume here that λ remains virtually constant during 2π/Ω.

2.1. Characteristic torque function

The thermal YORP torque by emitted photons is given by (NV07, NV08) $\mathit{\tau }\mathrm{=}\mathrm{-}\frac{\mathrm{2}\mathit{\gamma }}{\mathrm{3}{\mathit{v}}_{\mathit{c}}}{\mathrm{\int }}_{\mathrm{ℬ}}\mathrm{(}{x}\mathrm{\times }\mathit{\eta }\mathrm{)}{\mathit{T}}^{\mathrm{4}}\hspace{0.17em}\mathrm{d}\mathit{S,}$(5)where the integral is taken over the whole surface of the body ℬ (of arbitrary shape and topology), T is the surface temperature, γ is the product of the emissivity of the material in thermal wavelengths and the Stefan-Boltzmann constant, and v_c is the speed of light. (As shown in Nesvorný & Vokrouhlický 2008b, the torque due to absorbed photons vanishes identically when averaged over orbit revolution and body rotation.)

We define the characteristic torque function of a body as the three-dimensional vector $\mathrm{𝒯}\mathrm{\left(}\mathit{\xi }\mathrm{\right)}\mathrm{:}\mathrm{=}{x}\mathrm{\left(}\mathit{\xi }\mathrm{\right)}\mathrm{\times }\mathit{\eta }\mathrm{\left(}\mathit{\xi }\mathrm{\right)}\mathit{,}$(6)where ξ is any two-dimensional parametrization of the location on the surface. For a convex body, $\mathrm{d}\mathit{S}\mathrm{=}\mathit{G}\mathrm{\left(}\mathit{\eta }\mathrm{\right)}\hspace{0.17em}\mathrm{d}\mathit{\sigma ,} \mathrm{d}\mathit{\sigma }\mathrm{=}\sin \mathit{\vartheta }\hspace{0.17em}\mathrm{d}\mathit{\vartheta }\hspace{0.17em}\mathrm{d}\mathit{\psi ,}$(7)where G(η) is the curvature function, so we consider on S²: $\mathrm{𝒯}\mathrm{\left(}\mathit{\eta }\mathrm{\right)}\mathrm{:}\mathrm{=}\mathit{G}\mathrm{\left(}\mathit{\eta }\mathrm{\right)}{x}\mathrm{\left(}\mathit{\eta }\mathrm{\right)}\mathrm{\times }\mathit{\eta }$(8)(cf. BRV11). For asteroids, G(η) is typically an initial product from observational data, and x(η) is derived from it as the solution of the Minkowski problem (usually as a convex polyhedron; Kaasalainen et al. 1992, 2006).

Expanding as a Laplace series, we have $\mathrm{𝒯}\mathrm{\left(}\mathit{\eta }\mathrm{\right)}\mathrm{=}\sum _{\mathit{l}\hspace{0.17em}\mathrm{\ge }\hspace{0.17em}\mathrm{0}}\sum _{\mathit{m}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathrm{-}\mathit{l}}^{\mathit{l}}{{T}}_{\mathit{lm}}\hspace{0.17em}{\mathit{Y}}_{\mathit{l}}^{\mathit{m}}\mathrm{\left(}\mathit{\vartheta ,\psi }\mathrm{\right)}\mathit{.}$(9)Here we ignore the customary normalization constants of the spherical harmonics ${\mathit{Y}}_{\mathit{l}}^{\mathit{m}}$ as they can always be introduced as normalization factors when actually taking orthogonality integrals, so ${\mathit{Y}}_{\mathit{l}}^{\mathit{m}}\mathrm{\left(}\mathit{\vartheta ,\psi }\mathrm{\right)}\mathrm{=}{\mathit{P}}_{\mathit{l}}^{\mathit{m}}\mathrm{\left(}\cos \mathit{\vartheta }\mathrm{\right)}{\mathit{e}}^{\mathit{im\psi }}$. Written explicitly as a real-valued function, ${\mathrm{𝒯}}_{\mathit{i}}\mathrm{\left(}\mathit{\eta }\mathrm{\right)}\mathrm{=}\sum _{\mathit{l}\mathrm{\ge }\mathrm{0}}\sum _{\mathit{m}\mathrm{=}\mathrm{0}}^{\mathit{l}}{\mathit{P}}_{\mathit{l}}^{\mathit{m}}\mathrm{\left(}\mathit{\vartheta ,\psi }\mathrm{\right)}\mathrm{\left[}{\mathit{T}}_{\mathit{lm}}^{\mathit{i,}\cos }\cos \mathrm{\right(}\mathit{m\psi }\mathrm{)}\mathrm{+}{\mathit{T}}_{\mathit{lm}}^{\mathit{i,}\sin }\sin \mathrm{(}\mathit{m\psi }{\mathrm{\left)}}^{\mathrm{\right]}}\mathit{.}$(10)The possible symmetries of are evident in T_lm. For example, geometric inspection shows that if the body is symmetric about a line in the xy-plane, T_l0 ≡ 0. If the body is quadrant-symmetric (e.g., an ellipsoid), the Fourier terms are , with m odd, while with m even. The coefficients T_lm are readily calculated with the formulae in the appendix.

3. YORP torque and its components for K = 0

When there is no thermal inertia; i.e., the thermal conductivity K = 0, the absorbed flux is emitted immediately. Denoting $\mathit{\mu }\mathrm{:}\mathrm{=}\mathit{\eta }\mathrm{·}\mathit{\omega ,}$(11)the immediate emission means that, for a surface patch dS, $\mathit{\gamma }{\mathit{T}}^{\mathrm{4}}\mathrm{=}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{A}\mathrm{)}{\mathit{F}}_{\mathit{o}}\mathit{\mu ,}$(12)where A is the surface albedo and F_o the incoming solar flux. Thus we have, denoting the illuminated portion of the surface by A₊ and writing the convex version explicitly, $\mathit{\tau }\mathrm{=}\mathit{Q}{\mathrm{\int }}_{{\mathit{A}}_{\mathrm{+}}}\mathit{\mu }\hspace{0.17em}\mathrm{𝒯}\mathrm{\left(}\mathit{\xi }\mathrm{\right)}\hspace{0.17em}\mathit{dS}\mathrm{=}\mathit{Q}{\mathrm{\int }}_{\mathit{\mu }\hspace{0.17em}\mathrm{\ge }\hspace{0.17em}\mathrm{0}}\mathit{\mu }\hspace{0.17em}\mathrm{𝒯}\mathrm{\left(}\mathit{\eta }\mathrm{\right)}\hspace{0.17em}\mathrm{d}\mathit{\sigma ,}$(13)where Q has the physical dimension of pressure: $\mathit{Q}\mathrm{=}\mathrm{-}\frac{\mathrm{2}}{\mathrm{3}{\mathit{v}}_{\mathit{c}}}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{A}\mathrm{)}{\mathit{F}}_{\mathit{o}}\mathit{.}$(14)The effect of albedo variegation A(η) on the surface can be absorbed in by replacing G(η) with [1 − A(η)] G(η) in (8) and defining Q =  −2F₀/3v_c.

The YORP torque in the body frame is now obtained in the same way as the opposition brightness for convex bodies (Kaasalainen et al. 1992; now we just have the coefficients of a three-dimensional Laplace series). Thus we have the fundamental result for the YORP torque τ: $\mathit{\tau }\mathrm{\left(}\mathit{\theta ,\varphi }\mathrm{\right)}\mathrm{=}\mathit{Q}\sum _{\mathit{lm}}{{t}}_{\mathit{lm}}\hspace{0.17em}{\mathit{Y}}_{\mathit{l}}^{\mathit{m}}\mathrm{\left(}\mathit{\theta ,\varphi }\mathrm{\right)}\mathit{,}$(15)where, for convex bodies, t_lm = 2πk_l   T_lm, and ${\mathit{k}}_{\mathit{l}}\mathrm{=}\frac{\mathrm{1}}{\mathit{l}\mathrm{+}\mathrm{2}}\mathrm{(}\mathrm{-}\mathrm{1}{\mathrm{)}}^{\mathit{l}\mathit{/}\mathrm{2}\mathrm{-}\mathrm{1}}\frac{\mathrm{(}\mathit{l}\mathrm{-}\mathrm{3}\mathrm{)}\mathrm{!}\mathrm{!}}{\mathit{l}\mathrm{!}\mathrm{!}}\mathit{,} \mathit{l}\mathrm{\ge }\mathrm{2}\mathit{,}\hspace{0.17em}\hspace{0.17em}\mathit{l}\hspace{0.17em}\mathrm{even}\mathit{.}$(16)Other values of l are not of interest as the odd ones drop out for all shapes in averaging (see below), and T₀₀ vanishes identically due to Gauss’ theorem.

As pointed out in BRV11, the torque τ for any shape, computed on a grid of (θ,ϕ) by ray-tracing, can, of course, be expressed as a Laplace series by using, e.g., Lebedev-Laikov quadratures (Kaasalainen et al. 2012). The coefficients t_lm can be used in the following analytical motion-averaged formulae that no longer require convexity, although the relevant t_lm tend to be unstable, as discussed in Sect. 5.

3.1. Rotation speed

Averaging over one rotation, we have $\mathit{\tau }\mathrm{\left(}\mathit{ϵ,\lambda }\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{2}\mathit{\pi }}\mathit{\tau }\mathrm{\left(}\mathit{\theta ,\varphi }\mathrm{\right)}\hspace{0.17em}\mathrm{d}\mathit{\varphi }\mathrm{=}\mathit{Q}\sum _{\mathit{l}}{{t}}_{\mathit{l}\mathrm{0}}\hspace{0.17em}{\mathit{P}}_{\mathit{l}}\mathrm{\left(}\cos \mathit{\lambda }\sin \mathit{ϵ}\mathrm{\right)}\mathit{.}$(17)The secular torque is obtained by averaging over λ: $\mathit{\tau ̅}\mathrm{\left(}\mathit{ϵ}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{2}\mathit{\pi }}\mathit{\tau }\mathrm{\left(}\mathit{ϵ,\lambda }\mathrm{\right)}\hspace{0.17em}\mathrm{d}\mathit{\lambda }\mathrm{=}\mathit{Q}\sum _{\mathit{l}}{{t}}_{\mathit{l}\mathrm{0}}\hspace{0.17em}{\mathit{I}}_{\mathit{l}}\mathrm{\left(}\mathit{ϵ}\mathrm{\right)}\mathit{,}$(18)where ${\mathit{I}}_{\mathit{l}}\mathrm{\left(}\mathit{ϵ}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\mathit{\pi }}{\mathrm{\int }}_{\mathrm{0}}^{\mathit{\pi }}{\mathit{P}}_{\mathit{l}}\mathrm{\left(}\cos \mathit{\lambda }\sin \mathit{ϵ}\mathrm{\right)}\hspace{0.17em}\mathrm{d}\mathit{\lambda }\mathrm{=}\frac{\mathrm{1}}{\mathit{\pi }}{\mathrm{\int }}_{\mathrm{-}\sin \mathit{ϵ}}^{\sin \mathit{ϵ}}\frac{{\mathit{P}}_{\mathit{l}}\mathrm{\left(}\mathit{x}\mathrm{\right)}}{\sqrt{{\sin }^{\mathrm{2}}\mathit{ϵ}\mathrm{-}{\mathit{x}}^{\mathrm{2}}}}\hspace{0.17em}\mathrm{d}\mathit{x}\mathit{.}$(19)Thus I_l = 0 for odd l, and for even l we obtain, using the recursion relations of Legendre polynomials, ${\mathit{I}}_{\mathit{l}}^{\mathrm{even}}\mathrm{\left(}\mathit{ϵ}\mathrm{\right)}\mathrm{=}\frac{\mathrm{2}}{\mathit{\pi }}{\mathrm{\int }}_{\mathrm{0}}^{\sin \mathit{ϵ}}\frac{{\mathit{P}}_{\mathit{l}}\mathrm{\left(}\mathit{x}\mathrm{\right)}}{\sqrt{{\sin }^{\mathrm{2}}\mathit{ϵ}\mathrm{-}{\mathit{x}}^{\mathrm{2}}}}\hspace{0.17em}\mathrm{d}\mathit{x}\mathrm{=}{\mathit{C}}_{\mathit{l}}\hspace{0.17em}{\mathit{P}}_{\mathit{l}}\mathrm{\left(}\cos \mathit{ϵ}\mathrm{\right)}\mathit{,}$(20)where ${\mathit{C}}_{\mathit{l}}\mathrm{=}\mathrm{(}\mathrm{-}\mathrm{1}{\mathrm{)}}^{\mathit{l}\mathit{/}\mathrm{2}}\frac{\mathrm{(}\mathit{l}\mathrm{-}\mathrm{1}\mathrm{)}\mathrm{!}\mathrm{!}}{\mathit{l}\mathrm{!}\mathrm{!}}\mathrm{·}$(21)Our final expression for is now simply $\mathit{\tau ̅}\mathrm{\left(}\mathit{ϵ}\mathrm{\right)}\mathrm{=}\mathit{Q}\sum _{\mathit{l}\hspace{0.17em}\mathrm{even}}{\mathit{C}}_{\mathit{l}}\hspace{0.17em}{{t}}_{\mathit{l}\mathrm{0}}\hspace{0.17em}{\mathit{P}}_{\mathit{l}}\mathrm{\left(}\cos \mathit{ϵ}\mathrm{\right)}\mathit{,}$(22)which manifestly vanishes for azimuthally symmetric bodies due to the structure of T_lm discussed above. The torque component that alters the rotation speed of the body is the z-component given by the coefficients ${\mathit{t}}_{\mathit{l}\mathrm{0}}^{\mathrm{\left(}\mathrm{3}\mathrm{\right)}}$. If the series is dominated by ${\mathit{t}}_{\mathrm{20}}^{\mathrm{\left(}\mathrm{3}\mathrm{\right)}}$, this component essentially vanishes at the roots of P₂(cosϵ), i.e., at ϵ ≈ 55° and ϵ ≈ 125° (cf. NV07).

Shape-related phenomena that are apparent instabilities are obvious from this result. These are actually critical phenomena on a shape close to t_l0 = 0. Consider, say, a pointed local feature on an otherwise symmetric (e.g., ellipsoidal) body at ϑ = π/2, ψ = δ > 0. This feature produces a small . On the other hand, a similar arrangement symmetric about the x-axis, with the center of the feature at ψ =  −δ, yields exactly . Note that δ can be arbitrarily small. One can view this as a minuscule shift of a local feature on the body: practically invisible to the eye, yet the YORP torque changes its sign. This is an analogue of the numerical “boulder shift” phenomenon discussed in Statler (2009).

Fig. 1 Convex model shapes of asteroids Psyche (top) and Apollo (bottom).

Fig. 2 Secular YORP torque for the shapes of Fig. 1 respectively, computed analytically using Eq. (22) for values l ≤ 2, l ≤ 4 and l ≤ 6. Clearly, l = 2 is the dominant term, and the series converges rapidly for l > 2. The obliquity ϵ is given in radians.

As an example, we evaluated from Eq. (22) for the convex shapes shown in Fig. 1. In Fig. 2, the various linestyles for different numbers of ${\mathit{T}}_{\mathit{l}\mathrm{0}}^{\mathrm{\left(}\mathrm{3}\mathrm{\right)}}$-terms show the rapid convergence of the series. Note that, throughout this paper, we only analyze the shapes of asteroids, so we do not take sizes and masses into account (and set Q = 1). Thus the scale of is arbitrary.

3.2. Obliquity

The torque component τ_ϵ that alters the obliquity of the rotation axis is given by the rotational transformation of Eq. (4): ${\mathit{\tau }}_{\mathit{ϵ}}\mathrm{=}{\left[{{R}}_{\mathit{z}}\mathrm{(}\mathrm{-}{\mathit{\varphi }}^{\mathrm{\prime }}\mathrm{)}\mathit{\tau }\right]}_{\mathrm{1}}\mathrm{=}{\mathit{\tau }}_{\mathrm{1}}\cos {\mathit{\varphi }}^{\mathrm{\prime }}\mathrm{-}{\mathit{\tau }}_{\mathrm{2}}\sin {\mathit{\varphi }}^{\mathrm{\prime }}\mathit{.}$(23)In the same oblique coordinate frame (i.e., the z-axis is aligned with the rotation axis and the x-component of τ is τ_ϵ), the direction of the Sun is given by $\mathit{\omega }\mathrm{=}{{R}}_{\mathit{y}}\mathrm{\left(}\mathit{ϵ}\mathrm{\right)}\left(\begin{array}{c}\\ \cos \mathit{\lambda }\\ \sin \mathit{\lambda }\\ \mathrm{0}\end{array}\right)\mathrm{=}\left(\begin{array}{c}\\ \cos \mathit{ϵ}\cos \mathit{\lambda }\\ \sin \mathit{\lambda }\\ \sin \mathit{ϵ}\cos \mathit{\lambda }\end{array}\right)\mathrm{:}\mathrm{=}\left(\begin{array}{c}\\ \sin \mathit{\theta }\cos {\mathit{\varphi }}_{\mathit{s}}\\ \sin \mathit{\theta }\sin {\mathit{\varphi }}_{\mathit{s}}\\ \cos \mathit{\theta }\end{array}\right)\mathit{.}$(24)Thus ${\mathit{\varphi }}^{\mathrm{\prime }}\mathrm{=}{\mathit{\varphi }}_{\mathit{s}}\mathrm{-}\mathit{\varphi ,}$(25)and substituting this to (23) yields ${\mathit{\tau }}_{\mathit{ϵ}}\mathrm{=}\cos \mathit{\varphi }\mathrm{(}{\mathit{\tau }}_{\mathrm{1}}\cos {\mathit{\varphi }}_{\mathit{s}}\mathrm{-}{\mathit{\tau }}_{\mathrm{2}}\sin {\mathit{\varphi }}_{\mathit{s}}\mathrm{)}\mathrm{+}\sin \mathit{\varphi }\mathrm{(}{\mathit{\tau }}_{\mathrm{1}}\sin {\mathit{\varphi }}_{\mathit{s}}\mathrm{+}{\mathit{\tau }}_{\mathrm{2}}\cos {\mathit{\varphi }}_{\mathit{s}}\mathrm{)}\mathit{.}$(26)From Eq. (24) we find that $\cos {\mathit{\varphi }}_{\mathit{s}}\mathrm{=}\frac{\cos \mathit{ϵ}\cos \mathit{\lambda }}{\sqrt{\mathrm{1}\mathrm{-}{\sin }^{\mathrm{2}}\mathit{ϵ}{\cos }^{\mathrm{2}}\mathit{\lambda }}}\mathit{,}\sin {\mathit{\varphi }}_{\mathit{s}}\mathrm{=}\frac{\sin \mathit{\lambda }}{\sqrt{\mathrm{1}\mathrm{-}{\sin }^{\mathrm{2}}\mathit{ϵ}{\cos }^{\mathrm{2}}\mathit{\lambda }}}\mathit{,}$(27)so, substituting τ₁ and τ₂ from Eq. (15), using real-valued t_lm, and averaging over ϕ with Eq. (25), $\begin{array}{ccc}{\mathit{\tau }}_{\mathit{ϵ}}\mathrm{\left(}\mathit{\lambda ,ϵ}\mathrm{\right)}& \mathrm{=}& \\ & & \mathrm{+}\sin {\mathit{\varphi }}_{\mathit{s}}{\mathrm{(}{\mathit{t}}_{\mathit{l}\mathrm{1}}^{\mathrm{1}\mathit{,}\sin }\mathrm{-}{{\mathit{t}}_{\mathit{l}\mathrm{1}}^{\mathrm{2}\mathit{,}\cos }}^{\mathrm{)}}}^{\mathrm{]}}\mathit{,}\end{array}$(28)and averaging over λ, we obtain $\mathit{\tau ̅}\mathit{ϵ}\mathrm{\left(}\mathit{ϵ}\mathrm{\right)}\mathrm{=}\frac{\mathit{Q}}{\mathrm{2}}\sum _{\mathit{l}}\mathrm{(}{\mathit{t}}_{\mathit{l}\mathrm{1}}^{\mathrm{1}\mathit{,}\cos }\mathrm{+}{{\mathit{t}}_{\mathit{l}\mathrm{1}}^{\mathrm{2}\mathit{,}\sin }}^{\mathrm{)}}\hspace{0.17em}{\mathit{I}}_{\mathit{l}\mathrm{1}}\mathrm{\left(}\mathit{ϵ}\mathrm{\right)}\mathrm{+}\mathit{O,}$(29)where $\begin{array}{ccc}{\mathit{I}}_{\mathit{l}\mathrm{1}}\mathrm{\left(}\mathit{ϵ}\mathrm{\right)}& \mathrm{=}& \frac{\mathrm{1}}{\mathit{\pi }}{\mathrm{\int }}_{\mathrm{0}}^{\mathit{\pi }}{\mathit{P}}_{\mathit{l}}^{\mathrm{1}}\mathrm{\left(}\cos \mathit{\lambda }\sin \mathit{ϵ}\mathrm{\right)}\frac{\cos \mathit{ϵ}\cos \mathit{\lambda }}{\sqrt{\mathrm{1}\mathrm{-}{\sin }^{\mathrm{2}}\mathit{ϵ}{\cos }^{\mathrm{2}}\mathit{\lambda }}}\hspace{0.17em}\mathrm{d}\mathit{\lambda }\\ & \mathrm{=}& \frac{\mathrm{1}}{\mathit{\pi }}\frac{\cos \mathit{ϵ}}{\sin \mathit{ϵ}}{\mathrm{\int }}_{\mathrm{-}\sin \mathit{ϵ}}^{\sin \mathit{ϵ}}\frac{\mathit{x}\hspace{0.17em}{\mathit{P}}_{\mathit{l}}^{\mathrm{1}}\mathrm{\left(}\mathit{x}\mathrm{\right)}}{\sqrt{\mathrm{1}\mathrm{-}{\mathit{x}}^{\mathrm{2}}}}\frac{\mathrm{d}\mathit{x}}{\sqrt{{\sin }^{\mathrm{2}}\mathit{ϵ}\mathrm{-}{\mathit{x}}^{\mathrm{2}}}}\mathit{,}\end{array}$(30)so I_l1 = 0 for odd l, and O denotes the sinϕ_s-terms that vanish because of the parity of ${\mathit{P}}_{\mathit{l}}^{\mathrm{1}}$ and the requirement of an antisymmetric τ_ϵ: . Computing the integral in the same way as with , we have $\mathit{\tau ̅}\mathit{ϵ}\mathrm{\left(}\mathit{ϵ}\mathrm{\right)}\mathrm{=}\frac{\mathit{Q}}{\mathrm{2}}\sum _{\mathit{l}\hspace{0.17em}\mathrm{even}}{\mathit{C}}_{\mathit{l}}\hspace{0.17em}\mathrm{(}{\mathit{t}}_{\mathit{l}\mathrm{1}}^{\mathrm{1}\mathit{,}\cos }\mathrm{+}{{\mathit{t}}_{\mathit{l}\mathrm{1}}^{\mathrm{2}\mathit{,}\sin }}^{\mathrm{)}}\hspace{0.17em}{\mathit{P}}_{\mathit{l}}^{\mathrm{1}}\mathrm{\left(}\cos \mathit{ϵ}\mathrm{\right)}\mathit{.}$(31)Owing to the structure of T_lm, this vanishes for quadrant-symmetric bodies. In BRV11, Eqs. (29), (C14), and (C15) are alternative versions of Eqs. (15), (22), and (31) here.

4. Nonzero thermal conductivity

Let us define the insolation factor ins(μ,ξ) as ins(μ,ξ) = μ if the point ξ on the surface is illuminated, and otherwise ins(μ,ξ) = 0. For constant λ and ϵ, ins(μ,ξ) is obviously a cyclic function of ϕ′, so it can be expanded as a Fourier series: $\begin{array}{ccc}\mathrm{ins}\mathrm{\left(}\mathit{\mu ,\xi }\mathrm{\right)}& \mathrm{=}& \\ & \mathrm{=}& \sum _{\mathit{n}\mathrm{\ge }\mathrm{0}}{\mathit{a}}_{\mathit{n}}\mathrm{\left(}\mathit{\xi ,\lambda ,ϵ}\mathrm{\right)}\cos \mathit{n}{\mathit{\varphi }}^{\mathrm{\prime }}\mathrm{+}{\mathit{b}}_{\mathit{n}}\mathrm{\left(}\mathit{\xi ,\lambda ,ϵ}\mathrm{\right)}\sin \mathit{n}{\mathit{\varphi }}^{\mathrm{\prime }}\mathit{.}\end{array}$(32)The coefficients d_n(ξ,λ,ϵ) can always be readily computed, analytically for convex bodies and by ray-tracing and FFT for others. In this analysis, however, we do not need their actual values; it is sufficient to know that they exist and fully define ins(μ,ξ) (even when it is discontinuous due to shadowing by nonconvexities).

Writing τ by using ins(μ,ξ) instead of A₊, we have, for K = 0, $\begin{array}{ccc}\mathit{\tau }\mathrm{\left(}\mathit{\lambda ,ϵ,}{\mathit{\varphi }}^{\mathrm{\prime }}\mathrm{\right)}& \mathrm{=}& \\ & \mathrm{=}& \mathit{Q}\sum _{\mathit{n}}{\mathrm{e}}^{\mathrm{i}\mathit{n}{\mathit{\varphi }}^{\mathrm{\prime }}}{\mathrm{\int }}_{\mathrm{ℬ}}\mathrm{𝒯}\hspace{0.17em}\mathrm{\left(}\mathit{\xi }\mathrm{\right)}{\mathrm{d}}_{\mathit{n}}\mathrm{\left(}\mathit{\xi ,\lambda ,ϵ}\mathrm{\right)}\hspace{0.17em}\mathrm{d}\mathit{S}\mathit{.}\end{array}$(33)When averaging over ϕ′, the only contributing term is the one with n = 0.

For K ≠ 0, the heat diffusion equation $\mathit{\rho }{\mathit{c}}_{\mathit{p}}\frac{\mathit{\partial T}}{\mathit{\partial t}}\mathrm{=}\mathit{K}\frac{{\mathit{\partial }}^{\mathrm{2}}\mathit{T}}{\mathit{\partial }{\mathit{\xi }}^{\mathrm{2}}}\mathit{,}$(34)where ρ is the density of the material, c_p is the specific heat capacity, and the vertical direction ξ is aligned with η, can be solved with suitable physical approximations, boundary conditions, and a periodic ansatz of the form (32) (see NV08). With damping factors Ψ_n and phase lags Δφ_n defined by ${\mathrm{\Psi }}_{\mathit{n}}\mathrm{=}{\mathrm{(}\mathrm{1}\mathrm{+}\mathrm{2}{\mathrm{\Theta }}_{\mathit{n}}\mathrm{+}\mathrm{2}{{\mathrm{\Theta }}_{\mathit{n}}^{\mathrm{2}}}^{\mathrm{)}}}^{-1}\mathit{,} \mathrm{\Delta }{\mathit{\phi }}_{\mathit{n}}\mathrm{=}\mathrm{sgn}\mathrm{\left(}\mathit{n}\mathrm{\right)}\arctan \frac{{\mathrm{\Theta }}_{\mathit{n}}}{{\mathrm{\Theta }}_{\mathit{n}}\mathrm{+}\mathrm{1}}\mathit{,}$(35)where ${\mathrm{\Theta }}_{\mathit{n}}\mathrm{=}\frac{\sqrt{\mathit{\rho }{\mathit{c}}_{\mathit{p}}}}{\mathrm{4}\mathit{\gamma }{\mathit{T}}_{\mathrm{0}}^{\mathrm{3}}}\sqrt{\frac{\mathrm{1}}{\mathrm{2}}\mathrm{\left|}\mathit{n}\mathrm{\right|}\hspace{0.17em}\mathit{K}\hspace{0.17em}\mathrm{\Omega }}\mathit{,} {\mathit{T}}_{\mathrm{0}}^{\mathrm{4}}\mathrm{=}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{A}\mathrm{)}{\mathit{F}}_{\mathit{o}}\mathit{/}\mathit{\gamma ,}$(36)an approximate solution is ${\mathit{T}\mathrm{(}\mathit{t}{\mathrm{)}}^{\mathrm{4}}|}_{\mathit{\xi }\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathrm{0}}\mathrm{\simeq }\mathrm{(}\mathrm{1}\mathrm{-}\mathit{A}\mathrm{)}{\mathit{F}}_{\mathit{o}}{\mathit{\gamma }}^{-1}\sum _{\mathit{n}}{\mathrm{\Psi }}_{\mathit{n}}\hspace{0.17em}{\mathit{d}}_{\mathit{n}}\mathrm{\left(}\mathit{\eta ,\lambda ,ϵ}\mathrm{\right)}\hspace{0.17em}{\mathrm{e}}^{\mathrm{i}\mathrm{(}\mathit{n}{\mathit{\varphi }}^{\mathrm{\prime }}\mathrm{-}\mathrm{\Delta }{\mathit{\phi }}_{\mathit{n}}\mathrm{)}}\mathit{,}$(37)so $\mathit{\tau }\mathrm{\left(}\mathit{\lambda ,ϵ,}{\mathit{\varphi }}^{\mathrm{\prime }}\mathrm{\right)}\mathrm{=}\mathit{Q}\sum _{\mathit{n}}{\mathrm{\Psi }}_{\mathit{n}}\hspace{0.17em}{\mathrm{e}}^{\mathrm{i}\mathrm{(}\mathit{n}{\mathit{\varphi }}^{\mathrm{\prime }}\mathrm{-}\mathrm{\Delta }{\mathit{\phi }}_{\mathit{n}}\mathrm{)}}{\mathrm{\int }}_{\mathrm{ℬ}}\mathrm{𝒯}\mathrm{\left(}\mathit{\xi }\mathrm{\right)}{\mathit{d}}_{\mathit{n}}\mathrm{\left(}\mathit{\xi ,\lambda ,ϵ}\mathrm{\right)}\hspace{0.17em}\mathrm{d}\mathit{S}\mathit{.}$(38)Again, averaging over ϕ′ retains only the n = 0 -contribution, and Ψ₀ = 1, Δφ₀ = 0, so the result is exactly the same as with K = 0. This shows that K has no influence on the rotation-averaged τ.

4.1. Obliquity torque

Using the rotation-transform Eq. (23) for τ_ϵ together with the real-valued Fourier series in (32), we have, from Eqs. (23) and (33), $\begin{array}{ccc}{\mathit{\tau }}_{\mathit{ϵ}}\mathrm{\left(}\mathit{\lambda ,ϵ,}{\mathit{\varphi }}^{\mathrm{\prime }}\mathrm{\right)}& \mathrm{=}& \mathit{Q}\sum _{\mathit{n}}[\cos \mathit{n}{\mathit{\varphi }}^{\mathrm{\prime }}\cos {\mathit{\varphi }}^{\mathrm{\prime }}{\mathrm{\int }}_{\mathrm{ℬ}}{\mathrm{𝒯}}_{\mathrm{1}}\mathrm{\left(}\mathit{\xi }\mathrm{\right)}{\mathit{a}}_{\mathit{n}}\mathrm{\left(}\mathit{\xi ,\lambda ,ϵ}\mathrm{\right)}\hspace{0.17em}\mathrm{d}\mathit{S}\\ & & \mathrm{-}\hspace{0.17em}\sin \mathit{n}{\mathit{\varphi }}^{\mathrm{\prime }}\sin {\mathit{\varphi }}^{\mathrm{\prime }}{\mathrm{\int }}_{\mathrm{ℬ}}{\mathrm{𝒯}}_{\mathrm{2}}\mathrm{\left(}\mathit{\xi }\mathrm{\right)}{\mathit{b}}_{\mathit{n}}\mathrm{\left(}\mathit{\xi ,\lambda ,ϵ}\mathrm{\right)}\hspace{0.17em}\mathrm{d}\mathit{S}]\mathrm{+}\mathit{O,}\end{array}$(39)where O denotes cross-terms that vanish in averaging over ϕ′. Upon this averaging, we obtain an expression for K = 0 (Δφ_n = 0) in the Fourier-series representation: ${\mathit{\tau }}_{\mathit{ϵ}}\mathrm{\left(}\mathit{\lambda ,ϵ}\mathrm{\right)}\mathrm{=}\frac{\mathit{Q}}{\mathrm{2}}{\mathrm{\int }}_{\mathrm{ℬ}}\left[{\mathrm{𝒯}}_{\mathrm{1}}\mathrm{\left(}\mathit{\xi }\mathrm{\right)}{\mathit{a}}_{\mathrm{1}}\mathrm{\left(}\mathit{\xi ,\lambda ,ϵ}\mathrm{\right)}\mathrm{-}{\mathrm{𝒯}}_{\mathrm{2}}\mathrm{\left(}\mathit{\xi }\mathrm{\right)}{\mathit{b}}_{\mathrm{1}}\mathrm{\left(}\mathit{\xi ,\lambda ,ϵ}\mathrm{\right)}\right]\hspace{0.17em}\mathrm{d}\mathit{S,}$(40)which is thus equivalent to the ϕ-averaged τ_ϵ of Eq. (28), when assuming constant λ during the averaging.

Next we expand the above expression to K ≠ 0. With a real-valued Fourier series solution for T⁴, averaging τ_ϵ(λ,ϵ,ϕ′) from Eq. (38) over ϕ′ yields $\begin{array}{ccc}{\mathit{\tau }}_{\mathit{ϵ}}\mathrm{\left(}\mathit{\lambda ,ϵ}\mathrm{\right)}& \mathrm{=}& \frac{\mathit{Q}}{\mathrm{2}}{\mathrm{\Psi }}_{\mathrm{1}}[\cos \mathrm{\Delta }{\mathit{\phi }}_{\mathrm{1}}{\mathrm{\int }}_{\mathrm{ℬ}}\left({\mathrm{𝒯}}_{\mathrm{1}}{\mathit{a}}_{\mathrm{1}}\mathrm{-}{\mathrm{𝒯}}_{\mathrm{2}}{\mathit{b}}_{\mathrm{1}}\right)\hspace{0.17em}\mathrm{d}\mathit{S}\\ & & \mathrm{-}\hspace{0.17em}\sin \mathrm{\Delta }{\mathit{\phi }}_{\mathrm{1}}{\mathrm{\int }}_{\mathrm{ℬ}}\left({\mathrm{𝒯}}_{\mathrm{2}}{\mathit{a}}_{\mathrm{1}}\mathrm{+}{\mathrm{𝒯}}_{\mathrm{1}}{\mathit{b}}_{\mathrm{1}}\right)\hspace{0.17em}\mathrm{d}\mathit{S}]\mathit{.}\end{array}$(41)In other words, ${\mathit{\tau }}_{\mathit{ϵ}}\mathrm{\left(}\mathit{\lambda ,ϵ}\mathrm{\right)}\mathrm{=}{\mathrm{\Psi }}_{\mathrm{1}}\mathrm{[}\cos \mathrm{\Delta }{\mathit{\phi }}_{\mathrm{1}}{\mathit{\tau }}_{\mathit{ϵ}}\mathrm{\left(}\mathit{\tau }{\mathrm{|}}_{\mathrm{\Delta }\mathit{\phi }\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathrm{0}}\mathrm{\right)}\mathrm{+}\sin \mathrm{\Delta }{\mathit{\phi }}_{\mathrm{1}}{\mathit{\tau }}_{\mathit{ϵ}}{\mathrm{\left(}\mathit{\tau }{\mathrm{|}}_{\mathrm{\Delta }\mathit{\phi }\hspace{0.17em}\mathrm{=}\hspace{0.17em}\frac{\mathit{\pi }}{\mathrm{2}}}\mathrm{\right)}}^{\mathrm{]}}\mathit{,}$(42)where τ_ϵ(τ | _{Δφ    =    α}) denotes τ_ϵ(λ,ϵ) computed with a phase-lagged τ evaluated at ϕ′ − α instead of ϕ′; a phase lag of α = π/2 changes a₁ →  − b₁ and b₁ → a₁ in the Fourier series. Using ϕ → ϕ_s − (ϕ′ − π/2), i.e., ϕ → ϕ + π/2, in the fundamental Eq. (15) for τ when averaging over ϕ and λ, we obtain our final result: $\begin{array}{ccc}\mathit{\tau ̅}\mathit{ϵ}\mathrm{\left(}\mathit{ϵ}\mathrm{\right)}& \mathrm{=}& \\ & & \mathrm{-}\hspace{0.17em}\sin \mathrm{\Delta }{\mathit{\phi }}_{\mathrm{1}}{\mathrm{(}{\mathit{t}}_{\mathit{l}\mathrm{1}}^{\mathrm{1}\mathit{,}\sin }\mathrm{+}{{\mathit{t}}_{\mathit{l}\mathrm{1}}^{\mathrm{2}\mathit{,}\cos }}^{\mathrm{)}}}^{\mathrm{]}}\mathit{.}\end{array}$(43)This means that, at large values of KΩ, the YORP torque acquires a significant component that does not vanish even for ellipsoids, for example.

5. YORP stability

If we consider the coefficients t_lm in the basic torque Eq. (15), we can see that those determining the secular torques are precisely the ones that are usually vanishingly small compared to the momentary torque level. For example, an elongated ellipsoid produces large nonsecular torques, while the secular ones vanish. Thus even a small change in the shape of a body inevitably alters the secular torques with a much larger relative factor. If t_lm are determined numerically from the momentary torques with various illumination directions on S², the secular “trash coefficients” are easily affected by their errors, and can even be drowned in the numerical noise. This also helps to understand YORP instability as a fundamental and inevitable fact.

The analytical expressions in Sects. 3 and 4 provide a general set of convenient basis functions for the presentation of numerical YORP computations in a compact form. A robust approach for general shapes, especially those far from convex surfaces, is to compute the secular torques numerically with ray-tracing and averaging, and then fit series in P_l(cosϵ) or ${\mathit{P}}_{\mathit{l}}^{\mathrm{1}}\mathrm{\left(}\cos \mathit{ϵ}\mathrm{\right)}$, l even, to them. This one-dimensional procedure is more stable than the determination of the Laplace series of Eq. (15). One way to compute torques with K ≠ 0 is to evaluate the cases with emission phase lags of 0 and π/2 and then form their superposition by modulating with the phase lag Δφ₁ to conform with Eq. (43).

Many asteroids are well approximated by convex shapes (up to some level of resolution), so if two shapes differ only slightly, how different are their YORP properties? This question concerns not only the validity regime of an approximation, but the general stability properties of the YORP effect as well.

In this context, it is useful to normalize the YORP torques to be size-invariant, pertaining only to the shape of the object. Since the dimension of the torque integrand is that of volume, we define the dimensionless YORP capacity of a shape as the maximal average torque per unit volume and unit flux pressure (or torque when the volume V and the pressure are scaled to unity): $\mathrm{𝒴:}\mathrm{=}\frac{\max \mathrm{\left|}\mathit{\tau ̅}\mathrm{3}\mathrm{\right(}\mathit{ϵ}\mathrm{\left)}\mathrm{\right|}}{\mathit{QV}}\mathit{,} {\mathrm{𝒴}}_{\mathit{ϵ}}\mathrm{:}\mathrm{=}\frac{\max \mathrm{\left|}\mathit{\tau ̅}\mathit{ϵ}\mathrm{\right(}\mathit{ϵ}\mathrm{\left)}\mathrm{\right|}}{\mathit{QV}}\mathrm{·}$(44)The YORP capacity thus describes the amplitude of the -curve of the unity-scaled case. For most convex shapes, is fast to evaluate directly from ${\mathit{T}}_{\mathrm{20}}^{\mathrm{\left(}\mathrm{3}\mathrm{\right)}}$ with Eq. (22), and similarly with Eqs. (31) or (43). Conversely, given for any shape, one can quickly compute the actual maximal torque by multiplying by QV.

To establish a scale for , consider a “propeller” shape consisting of two identical vertically symmetric wedges attached to each other inversion-symmetrically (Fig. 3). The maximal YORP torque is attained at ϵ = 0,π. The cross-section of a wedge is an isosceles triangle with height h and width w, and the length of a wedge is l. If we assume a long and thin body with h/l ≪ 1 (and w/l ≪ 1), we can neglect the torques from the ends of the wedges and terms of order (h/l)² and higher, so $\mathit{\tau ̂}\mathrm{3}\mathrm{\left(}\mathrm{0}\mathrm{\right)}\mathrm{\approx }\mathrm{𝒲}\frac{\mathrm{2}}{\mathit{\pi }}{\mathrm{\int }}_{\mathrm{0}}^{\mathit{\pi }\mathit{/}\mathrm{2}}\mathit{hl}\frac{\mathit{l}}{\mathrm{2}}\cos \mathit{\phi }\hspace{0.17em}\mathit{d\phi }\mathrm{=}\mathrm{𝒲}\frac{\mathit{h}{\mathit{l}}^{\mathrm{2}}}{\mathit{\pi }}\mathit{,}$(45)where the efficiency factor of the sharpness of the wedge is $\mathrm{𝒲}\mathrm{=}\mathrm{1}\mathrm{-}\frac{\mathrm{1}}{\sqrt{\mathrm{1}\mathrm{+}\mathrm{(}\mathrm{2}\mathit{w}\mathit{/}\mathit{h}{\mathrm{)}}^{\mathrm{2}}}}\mathrm{·}$(46)Thus, with V = hwl and estimating from , $\mathrm{𝒴}\mathrm{\approx }\frac{\mathrm{𝒲}}{\mathit{\pi }}\frac{\mathit{l}}{\mathit{w}}\mathrm{\approx }\mathrm{0.1}\frac{\mathit{l}}{\mathit{h}}$(47)for a flattening ratio w/h = 3/2 realistic for an asteroid. An infinitely long and thin body yields , so there is no upper bound for the YORP capacity in principle. For realistic shapes l/h ~ 1, we can estimate an upper bound .

Fig. 3 Wedge-propeller shape strongly prone to the rotational YORP torque, shown as a solid object and as a wire frame.

Any suitably asymmetric, infinitely long and thin shape yields , so there is no unique answer to the question “what is the most YORP-prone shape”, nor is there a unique solution to the problem “what is the most YORP-efficient realistic asteroid shape” since the term “realistic” can be defined in a number of ways. The wedge-propeller shape is more or less the optimal one for YORP-efficiency as it portrays the best “windmill asymmetry” profile.

To study the stability problem, let us estimate an upper bound for the change due to the perturbation Δx that affects the characteristic torque function: $\mathrm{𝒯}\mathrm{+}\mathrm{\Delta 𝒯}\mathrm{=}\mathrm{(}\mathit{G}\mathrm{+}\mathrm{\Delta }\mathit{G}\mathrm{)}\mathrm{(}{x}\mathrm{+}\mathrm{\Delta }{x}\mathrm{)}\mathrm{\times }\mathit{\eta }\mathit{.}$(48)Consider, e.g., the perturbed unit sphere, so that (the estimate is essentially similar for all convex bodies, but we choose the sphere for a simple example). Thus $\mathrm{\Delta }{{T}}_{\mathit{lm}}\mathrm{\approx }{\mathit{N}}_{\mathit{lm}}^{-1}{\mathrm{\int }}_{{\mathit{S}}^{\mathrm{2}}}\mathrm{\Delta }{x}\mathrm{\times }\mathit{\eta }{\mathit{Y}}_{\mathit{l}}^{\mathrm{\ast }\mathit{m}}\mathrm{\left(}\mathit{\eta }\mathrm{\right)}\hspace{0.17em}\mathrm{d}\mathit{\sigma ,}$(49)and we can approximate, with a perturbation Δx = f(η), where f(η) is a Laplace series such that x + Δx is convex and ⟨  ∥ f ∥  ⟩ _η = δ, $\mathrm{\left|}\mathrm{\Delta 𝒴}\mathrm{\right|}\mathrm{~}\mathrm{2}\mathit{\pi }\hspace{0.17em}\left|{\mathit{k}}_{\mathrm{2}}{\mathit{C}}_{\mathrm{2}}\hspace{0.17em}{\mathit{N}}_{\mathrm{20}}^{-1}\frac{\mathrm{3}}{\mathrm{4}\mathit{\pi }}{\mathrm{\int }}_{{\mathit{S}}^{\mathrm{2}}}{\left(\mathrm{\Delta }{x}\mathrm{\times }\mathit{\eta }\right)}_{\mathrm{3}}{\mathit{P}}_{\mathrm{2}}\mathrm{\left(}\mathit{\vartheta }\mathrm{\right)}\hspace{0.17em}\mathrm{d}\mathit{\sigma }\right|\lesssim \mathit{\delta }\mathit{/}\mathrm{10}\mathit{,}$(50)where the upper bound for the integral is estimated by applying the Legendre recursion formulae to the products with the components of η and using the orthogonality integral. A perturbation of δ = 0.1 can thus produce a change roughly a tenth of the wedge propeller’s above; in practice, the change is usually smaller than this.

When estimating the change of YORP properties due to perturbations, we only consider the change in amplitude. Indeed, taking the leading Legendre term as in Eq. (50) implies that the changes in the phase (the value of ϵ at the maximum) and the shape of are, for small shape perturbations, negligible compared to the amplitude change.

In addition to by surface fluctuations of a given size (that typically describe the convexity/nonconvexity or resolution difference between two shape models), we can consider a couple of other mechanisms for comparison: i) the “boulder effect” and ii) the effect of scale-invariant surface roughness. If we attach the wedge above to the surface of a sphere of radius r, and fix h = w = l = δr, we have , so adding a boulder of size r/10 changes by . If we consider a special arrangement of infinitesimal wedges on the surface, all pointing in an ordered way such that l/h = c and h² → dσ and hl → c   dσ, we obtain an upper bound $\frac{\mathrm{3}}{\mathrm{4}\mathit{\pi }{\mathit{r}}^{\mathrm{3}}}\frac{\mathrm{1}}{\mathit{\pi }}{\mathrm{\int }}_{\mathrm{0}}^{\mathit{\pi }}{\mathrm{\int }}_{\mathrm{0}}^{\mathit{\pi }\mathit{/}\mathrm{2}}\mathit{r}\sin \mathit{\theta }\mathrm{𝒲}\cos \mathit{\phi }{\mathit{r}}^{\mathrm{2}}\sin \mathit{\theta }\hspace{0.17em}\mathrm{d}\mathit{\theta }\hspace{0.17em}\mathrm{d}\mathit{\phi }\mathrm{=}\frac{\mathrm{3}\mathit{c}\mathrm{𝒲}}{\mathrm{8}\mathit{\pi }}\mathit{,}$(51)so a moderate surface roughness c ≈ 0.1 yields . In principle, one can coat a sphere such that it looks smooth and yet is, paradoxically, quite prone to YORP. The upper bound cannot be reached by a realistic unordered surface structure, but is plausible. Thus, taking a single large boulder or the unseen small-scale surface structure into account on a small asteroid can change up to .

Fig. 4 Nonconvex perturbation based on the model shape of the asteroid Psyche in Fig. 1. Here the value of the YORP capacity of the original body is 1.84 × 10-3, whereas the perturbed shape has . The relative radius perturbation ⟨ Δr ⟩ / ⟨ r ⟩ is 0.040.

Fig. 5 Plot of ⟨ Δr ⟩ / ⟨ r ⟩ versus for different global perturbations of the asteroid Apollo. A bounding line is drawn with the slope kmax = 0.22. As the relative perturbation ⟨ Δr ⟩ / ⟨ r ⟩ grows, the absolute change in remains bounded and well below the linear upper bound estimate. An alternative bounding line with 80% of the values of remaining under the line is plotted in the same figure, with the slope k80% = 0.13.

Given the analytical scale estimates of and above, are there YORP-stable asteroid shapes in reality? We measured the YORP values and their perturbations for some targets, when the shape perturbations were global nonconvex radius fluctuations Δr of varying magnitude (by low-order Laplace series), and was computed by ray-tracing. A sample perturbation is shown in Fig. 4. For each shape and sequence of perturbations , one can define an upper bound k_max such that $\mathrm{\left|}\mathrm{\Delta 𝒴}\mathrm{\right|}\lesssim {\mathit{k}}_{\max }\frac{\mathrm{⟨}\mathrm{\Delta }\mathit{r}\mathrm{⟩}}{\mathrm{⟨}\mathit{r}\mathrm{⟩}}\mathrm{·}$(52)To estimate this, we plotted ⟨ Δr ⟩ / ⟨ r ⟩ vs. for each shape; a sample case is shown in Fig. 5. A bounding line is approximately at k_max ≈  [0.1,0.2], consistent with the upper bound δ/10 estimated above. What is more, the slope k_max decreases as ⟨ Δr ⟩ / ⟨ r ⟩ grows ( may even saturate) so the linear upper bound estimate is only a rough approximation for small perturbations. For most shapes, k_max ~ 0.1 for a practical estimate. The estimate is not strict, nor can a shape perturbation be properly described by one figure ⟨ Δr ⟩ / ⟨ r ⟩ only.

The global perturbation scheme depicts the case of two alternative representations of the same body – for example, a lightcurve-based model and a radar-based one, or two versions of a model based on flyby data. If the model error ⟨ Δr ⟩ / ⟨ r ⟩ remains small, we can estimate shapes with (such as Apollo) to be reasonably YORP-stable against global model errors in the sense of the relative YORP error . Any instability in this error is caused by the inevitable smallness of rather than the size of .

Another interesting case is the “boulder effect”; it represents a local model error, or an actual local change of the shape due to some physical perturbation. We made random local changes to bodies and measured the YORP change; a scatter plot similar to that of global perturbations is shown in Fig. 6. In this case, most of the YORP changes are negligible, while some jump to . As in the analytical estimate above, the relative YORP changes remain well bounded for shapes with (such as the plotted case of Apollo), while shapes with are no longer YORP-stable against local shape perturbations.

In Table 1, we show sample values of for a number of asteroid shapes (obtained from the DAMIT website). Note again that is size-independent, and the samples represent typical asteroid shapes rather than estimates for the actual bodies, many of which are large and not affected by YORP.

Fig. 6 Plot of ⟨ Δr ⟩ / ⟨ r ⟩ versus for different local perturbations of the asteroid Apollo. Only one facet group of the asteroid is perturbed at a time, and only this section defines ⟨ Δr ⟩ / ⟨ r ⟩. With local perturbations less than 10%, the relative change in mostly remains below 10%. Therefore, we can conclude that the asteroid Apollo is locally YORP-stable.

To sum up the analytical and numerical stability estimates above as a rough order-of-magnitude classification scheme, we can define YORP stability classes as $\mathrm{stability}\hspace{0.17em}\hspace{0.17em}\mathrm{class}\mathrm{=}\mathit{n}\mathrm{=}\mathrm{round}\left({\log }_{\mathrm{10}}\mathrm{𝒴}\mathrm{+}\mathrm{1.8}\right)\mathit{,}$(53)where round(x) gives the integer closest to x, and the constant 1.8 was chosen to set the lower class limit to . Most asteroid shapes are in the regime n ≤  − 1; some reach n = 0, and mostly propeller-type shapes have n ≥ 1. Depending on the stability class, all the perturbation scenarios above can change from tens of percent to orders of magnitude for typical asteroid shapes in a quite natural manner. As a ballpark estimate, we could call shapes with n = 0 locally YORP-stable and globally semistable, while those with n =  −1 are locally semistable and globally unstable. All shapes in the classes n <  − 1 are YORP-unstable, and those in n ≥ 1 YORP-stable.

Most of the shapes in Table 1 are in n =  −1; only “Massalia”, “Apollo”, “Geographos”, and “Eger” are in n = 0. For the last three of these more stable ones, the YORP effect has been observed directly from lightcurves (Kaasalainen et al. 2007; Ďurech et al. 2008, 2012). The stability estimate is consistent with the good correspondence between the computed and observed YORP acceleration for Apollo and Geographos. For Eger, the correspondence is not as good. The value of thus indicates that this is more due to its unknown sizable nonconvex features than the smaller-scale irregularities (the cited value is based on the adopted nonconvex shape solution; the convex one has an even higher ). Indeed, Eger is notable for being one of the few asteroids for which lightcurves indicate strong nonconvexities (Ďurech & Kaasalainen 2003). Some main features can be sketched from photometric data (Ďurech et al. 2012), but the solution is not unique, which may cause the discrepancy between the expected and real YORP strengths.

As in the analytical estimates, the amplitude change dominated our numerical perturbation cases. This means that, for most of the -curve, determines the stability class. However, when ϵ is close to the values at which vanishes (i.e., usually around ϵ ≈ 55° and ϵ ≈ 125°), the changes in the phase and shape of the curve induced by the shape perturbation are significant for the relative change in . Near these obliquities, all shapes are YORP-unstable.

6. Conclusions and discussion

We have presented a concise formalism for analyzing YORP-related quantities. Using the characteristic torque function of the body, defined by Eq. (8), one can readily write the general Eq. (15) for the YORP torque as well as its main secular components of interest: the one changing the rate of rotation and the one altering the obliquity. The motion-averaged results can be given for both zero (Eqs. (22), (31)) and nonzero (Eq. (43)) thermal inertia. Since the averaging integrals are analytical regardless of the shape, the formulae also provide a convenient set of basis functions for representing YORP computations of nonconvex shapes.

The YORP effect as such is stable in principle; i.e., a small change in the shape leads to a small change in the YORP capacity . What causes the unpredictability is the relative change because most asteroids tend to have close to zero. Here we have derived a number of analytical and numerical estimates on YORP stability, and, based on these, surmise that most asteroids tend to be YORP-unstable, though some can reach a stable class. The stability class is determined by the YORP capacity : the larger the capacity, the better the relative YORP stability. For convex polyhedra, a simple approximation of is obtained from the P₂-coefficient of : $\mathrm{𝒴}\mathrm{\approx }\frac{\mathrm{5}}{\mathrm{64}\mathit{V}}|\sum _{\mathit{i}}{\mathit{A}}_{\mathit{i}}{\left({{x}}_{\mathit{i}}\mathrm{\times }{\mathit{\eta }}_{\mathit{i}}\right)}_{\mathrm{3}}\mathrm{(}\mathrm{3}{\cos }^{\mathrm{2}}{\mathit{\vartheta }}_{\mathit{i}}\mathrm{-}{\mathrm{1}}^{\mathrm{)}}|\mathit{,}$(54)where A_i, x_i, and η_i are the area, centroid, and unit normal of a facet, and ϑ_i is the polar angle of the normal (see Appendix). This approximation can also be used directly for slightly nonconvex shapes to obtain a first estimate; for these, the sum over facets is almost similar to a convex case.

Here we have simplified the YORP business to the extreme at the expense of accuracy, but Eq. (54) provides some equational economy: with it (and its variants for other torques), we can approximate and and estimate the YORP stability of a body. For more strongly nonconvex shapes, one can construct the best-fit convex shape by, e.g., scaling the convex hull suitably (semianalytical generalizations of NV07 via a perturbation of this shape are also relatively straightforward, but they lose the conceptual simplicity with no certainty of better accuracy). If the fluctuations around this surface do not exceed, say, a tenth part of the radius, the analytical result is reasonably reliable if . For , the accuracy of any YORP estimate (not just the types discussed here) is uncertain due to the instability. For , all YORP estimates are likely to be invalid, and only upper bounds can be given. The same applies to YORP estimates made near the obliquity at which the averaged net torque vanishes, regardless of .

An unstable class does not mean that all shape perturbations induce a large change in , but that suitable ones can do so. It would be interesting to explore the YORP stability properties in “shape space” by evaluating and for a large number of both convex and nonconvex shapes and their perturbations by craters and valleys. This would shed more light on the role of the YORP effect in the evolutionary history of asteroids. For example, a small asteroid with a low YORP capacity is likely to undergo a more stochastic YORP-affected spin evolution than one with a higher value of , since even small alterations of the shape may be significant for the YORP effect.

Acknowledgments

This work was supported by the Academy of Finland project “Modelling and applications of stochastic and regular surfaces in inverse problems”. We thank David Vokrouhlický for valuable comments that improved the manuscript.

References

Appendix A: Polyhedral representations

The torque over a planar facet of area A_i is given by the torque density at the centroid of the facet multiplied by the facet area, as is easily verified by taking the integral (13) on the facet alone in a coordinate system transformed to the facet plane. Thus, e.g., the YORP torque of a convex polyhedron at K = 0 is $\mathit{\tau }\mathrm{=}\mathit{Q}\sum _{\mathit{i}}\frac{{\mathit{A}}_{\mathit{i}}}{\mathrm{3}}{\mathit{\mu }}_{\mathit{i}}\sum _{\mathit{j}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathrm{1}}^{\mathrm{3}}{{x}}^{\mathrm{\left(}\mathit{ij}\mathrm{\right)}}\mathrm{\times }{\mathit{\eta }}_{\mathit{i}}\mathit{,}$(A.1)where x^(ij) are the vertices of the triangle facets, and their areas and outward unit normals are A_i and η_i.

Consider next ${{T}}_{\mathit{lm}}\mathrm{=}{\mathit{N}}_{\mathit{lm}}^{-1}{\mathrm{\int }}_{{\mathit{S}}^{\mathrm{2}}}\mathit{G}\mathrm{\left(}\mathit{\eta }\mathrm{\right)}\mathrm{(}{x}\mathrm{\times }\mathit{\eta }\mathrm{)}{\mathit{Y}}_{\mathit{l}}^{\mathrm{\ast }\mathit{m}}\mathrm{\left(}\mathit{\eta }\mathrm{\right)}\hspace{0.17em}\mathrm{d}\mathit{\sigma ,}$(A.2)where ${\mathit{N}}_{\mathit{lm}}\mathrm{=}\frac{\mathrm{4}\mathit{\pi }}{\mathrm{2}\mathit{l}\mathrm{+}\mathrm{1}}\frac{\mathrm{(}\mathit{l}\mathrm{+}\mathit{m}\mathrm{)}\mathrm{!}}{\mathrm{(}\mathit{l}\mathrm{-}\mathit{m}\mathrm{)}\mathrm{!}}$(A.3)for the unnormalized ${\mathit{Y}}_{\mathit{l}}^{\mathit{m}}\mathrm{\left(}\mathit{\vartheta ,\psi }\mathrm{\right)}\mathrm{=}{\mathit{P}}_{\mathit{l}}^{\mathit{m}}\mathrm{\left(}\cos \mathit{\vartheta }\mathrm{\right)}{\mathrm{e}}^{\mathrm{i}\mathit{m\psi }}$. In the real-valued version of the spherical harmonics series, the coefficients for the sine and cosine terms with m > 0 are twice the ones in the complex-valued case. The torque from Eq. (A.1) can be

interpreted as the Dirac-delta limit of the function-based form of Eq. (13), so we obtain an exact result for the coefficients T_lm of a convex polyhedron by applying the same limit to Eq. (A.2): ${{T}}_{\mathit{lm}}\mathrm{=}{\mathit{N}}_{\mathit{lm}}^{-1}\sum _{\mathit{i}}{\mathit{A}}_{\mathit{i}}\mathrm{(}{{x}}_{\mathit{i}}\mathrm{\times }{\mathit{\eta }}_{\mathit{i}}\mathrm{)}{\mathit{Y}}_{\mathit{l}}^{\mathrm{\ast }\mathit{m}}\mathrm{\left(}{\mathit{\eta }}_{\mathit{i}}\mathrm{\right)}\mathit{,}$(A.4)where x_i denotes the centroid of the facet i. If G(η) and x(η) are defined with spherical harmonics series, τ or T_lm can be computed analytically as mentioned earlier, or with quadratures, especially Lebedev-Laikov on S² (Kaasalainen et al. 2012).

The change in due to shifting the origin to some x₀; i.e., x → x − x₀, is easily computed by writing with $\mathrm{\Delta 𝒯}\mathrm{\left(}\mathit{\eta }\mathrm{\right)}\mathrm{=}\mathit{G}\mathrm{\left(}\mathit{\eta }\mathrm{\right)}\hspace{0.17em}{{x}}_{\mathrm{0}}\mathrm{\times }\mathit{\eta }\mathit{.}$(A.5)Thus the coefficients T_lm change as T_lm → T_lm − ΔT_lm, with $\begin{array}{ccc}\mathrm{\Delta }{{T}}_{\mathit{lm}}& \mathrm{=}& {\mathit{N}}_{\mathit{lm}}^{-1}{\mathrm{\int }}_{{\mathit{S}}^{\mathrm{2}}}\mathit{G}\mathrm{\left(}\mathit{\eta }\mathrm{\right)}\hspace{0.17em}\mathrm{(}{{x}}_{\mathrm{0}}\mathrm{\times }\mathit{\eta }\mathrm{)}{\mathit{Y}}_{\mathit{l}}^{\mathrm{\ast }\mathit{m}}\mathrm{\left(}\mathit{\eta }\mathrm{\right)}\hspace{0.17em}\mathrm{d}\mathit{\sigma }\\ & \mathrm{=}& {\mathit{N}}_{\mathit{lm}}^{-1}\sum _{\mathit{i}}{\mathit{A}}_{\mathit{i}}\mathrm{(}{{x}}_{\mathrm{0}}\mathrm{\times }{\mathit{\eta }}_{\mathit{i}}\mathrm{)}{\mathit{Y}}_{\mathit{l}}^{\mathrm{\ast }\mathit{m}}\mathrm{\left(}{\mathit{\eta }}_{\mathit{i}}\mathrm{\right)}\mathit{,}\end{array}$(A.6)which, if G(η) is written as a spherical harmonics series, is simple to evaluate by writing the products of ${\mathit{Y}}_{\mathit{l}}^{\mathit{m}}$ and the components of η with the Legendre recursion relations and using the orthogonality integral.

All Tables

All Figures

	Fig. 1 Convex model shapes of asteroids Psyche (top) and Apollo (bottom).
In the text

	Fig. 2 Secular YORP torque for the shapes of Fig. 1 respectively, computed analytically using Eq. (22) for values l ≤ 2, l ≤ 4 and l ≤ 6. Clearly, l = 2 is the dominant term, and the series converges rapidly for l > 2. The obliquity ϵ is given in radians.
In the text

	Fig. 3 Wedge-propeller shape strongly prone to the rotational YORP torque, shown as a solid object and as a wire frame.
In the text

	Fig. 4 Nonconvex perturbation based on the model shape of the asteroid Psyche in Fig. 1. Here the value of the YORP capacity of the original body is 1.84 × 10-3, whereas the perturbed shape has . The relative radius perturbation ⟨ Δr ⟩ / ⟨ r ⟩ is 0.040.
In the text

Fig. 5 Plot of ⟨ Δr ⟩ / ⟨ r ⟩ versus for different global perturbations of the asteroid Apollo. A bounding line is drawn with the slope kmax = 0.22. As the relative perturbation ⟨ Δr ⟩ / ⟨ r ⟩ grows, the absolute change in remains bounded and well below the linear upper bound estimate. An alternative bounding line with 80% of the values of remaining under the line is plotted in the same figure, with the slope k80% = 0.13.

In the text

	Fig. 6 Plot of ⟨ Δr ⟩ / ⟨ r ⟩ versus for different local perturbations of the asteroid Apollo. Only one facet group of the asteroid is perturbed at a time, and only this section defines ⟨ Δr ⟩ / ⟨ r ⟩. With local perturbations less than 10%, the relative change in mostly remains below 10%. Therefore, we can conclude that the asteroid Apollo is locally YORP-stable.
In the text