Analytic model for tangential YORP
Oleksiy Golubov


Oleksiy Golubov
oleksiy.golubov@karazin.ua
0000-0002-5720-2282]Oleksiy Golubov
Department of Aerospace Engineering Sciences, University of Colorado at Boulder 

429 UCB, Boulder, CO, 80309, USA
Karazin Kharkiv National University, 4 Svobody Sq., Kharkiv, 61022, Ukraine
Institute of Astronomy of V. N. Karazin Kharkiv National University, 35 Sumska Str., Kharkiv, 61022, Ukraine


The tangential YORP effect (TYORP) plays a significant role in the dynamical evolution of asteroids, and up to now has only been studied numerically.
This paper describes the first analytic model of the TYORP effect.
Although the model rests on numerous physical and mathematical simplifications,
the final analytic expression for TYORP is found to be in agreement 
with the results of rigorous numeric simulations to the accuracy of several tens per cent.

The obtained analytic expression is used to estimate the TYORP produced by the non-flat surface of regolith, – a contribution to TYORP that has never been considered. It is found that the contribution to TYORP arising from regolith can be comparable to the conventional TYORP produced by boulders.

Then, the analytic expression is fitted with a log-normal function
and used to integrate TYORP over all boulder sizes.
The general trend of TYORP for multiple boulders appears qualitatively similar to the trend of one boulder,
and also demonstrates a maximal TYORP at some particular rotation rate.
The obtained expression for integrated TYORP may be instrumental for simulations of evolution of asteroids subject to TYORP.

To conclude, the physical origin of TYORP is discussed in light of the constructed analytic model.







§ INTRODUCTION

The tangential YORP effect, or TYORP, appears when stones on the surface of an asteroid emit different amounts of infrared light eastward and westward, thus experiencing a net recoil force tangential to the asteroid's surface.
Until now, this effect has only been studied in numeric simulations <cit.>.
Although it was generally understood that the effect was due to the non-linearity of the heat emission law,
the detailed physics of the effect remained obscure.
Moreover, the question remained whether the whole effect could be attributed to numeric artifacts.

In this article I propose a minimalistic analytic model of the effect, which is based on the following simplifying assumptions:

1. Instead of solving partial differential equation for the heat conduction in a boulder, the boulder is split into two parts,
the eastern part and the western part, and the mean temperature of each part is introduced. Then the mathematical description of the model boils down to a system of two ordinary differential equations for the mean temperatures.

2. The incoming solar energy as a function of time is approximated by the sum of its zeroth and first order Fourier terms; 
all higher-order terms are neglected.

3. The first order Fourier term is treated perturbatively, as if it were small compared to the zeroth order term.

After these simplifications the problem can be easily solved analytically,
and the TYORP drag can be calculated.
The result fits the numeric simulations surprisingly well.

In Section <ref>, I provide a derivation of the analytic expression for TYORP.
In Section <ref>, the derived analytic expression is applied to different geometries, to test it and to make some new predictions.
In Section <ref>, I simplify the analytic expression and integrate it over different boulder sizes, to evaluate the total TYORP experienced by an asteroid.
In Section <ref>, I discuss how the derived analytic expression helps to better understand the physics of TYORP.	



§ GENERAL THEORY



 §.§ Derivation of the heat conduction equations

Heat balance within any volume part V of a boulder is governed by the following heat conduction equation in the integral form,

    Cρ∫_V∂ T/∂ tdV=
     κ∫_S_st∂ T/∂ X_idS_i+ 
     
     +(1-A)∫_S_spI_idS_i-
     ϵσ∫_S_spT^4dS .

Here T stands for the temperature.
The left-hand side describes the total heat energy increase in the volume V,
while the right-hand side is the sum of the heat conduction into this volume,
the direct solar heat absorbed by its open surface, 
and the negative heat emitted by the open surface.
The surface areas S_st, S_sp, and S_reg are parts of the volume's boundary bordering respectively
stone, space, and regolith, so that S_st+S_sp+S_reg=∂ V is the full boundary of the volume V 
(see the left-hand panel of Figure <ref>).
The heat conductivity of the stone is κ, 
its heat capacity is C, the density is ρ, 
the hemispherical albedo is A, and the emissivity is ϵ.
The heat conductivity of the regolith is assumed to be zero.
σ is Stefan–Boltzmann's constant, and 𝐈 is the vector of the incoming solar energy flux.



Now I nondimensionalize the variables.
Instead of time t, I use the rotation phase ϕ=ω t, with ω being the angular velocity of the asteroid.
By definition, ϕ=0 at noon.
The characteristic scales of length and temperature are

    L_cond=κ/((1-A)Φ)^3/4(ϵσ)^1/4 ,


    T_0=√((1-A)Φ/ϵσ) ,

with Φ being the solar constant.
Here T_0 is the equilibrium temperature at the subsolar point, while L_cond is the distance at which the temperature difference T_0 creates heat flux equal to AΦ.
I use these two scales to introduce the dimensionless variables x_i=X_i/L_cond and τ=T/T_0.
The relative importance of heat conduction with respect to heat absorption and emission is characterized by the thermal parameter

    θ=(Cρκω)^1/2/((1-A)Φ)^3/4(ϵσ)^1/4.


With these definitions, Eqn. (<ref>) transforms into

    θ^2∫_v∂τ/∂ϕdv=
     ∫_s_st∂τ/∂ x_ids_i+
     ∫_s_spI_i/Φds_i-
     ∫_s_spτ^4ds .

Here v and s correspond to the same volumes and areas as before,
but measured in the dimensionless variables x_i
instead of the dimensional variables X_i.

Now I separate the boulder into the western and the eastern parts and apply Eqn. (<ref>) to each part separately.
I assume that the boulder is symmetric, with the western and the eastern parts being mirror reflections of each other.
Let v henceforth denote the dimensionless volume of each half of the boulder (either v_w or v_e, see the right-hand panel of Figure <ref>),
and s_sp denote their equal dimensionless surface areas bordering space (either s_sp e or s_sp w in the right-hand panel of Figure <ref>).
I denote the mean dimensionless temperatures of the western and the eastern parts of the boulder via τ_w and τ_e correspondingly,

    τ_w   =   1/v∫_v_wτ  dv  , 
    τ_e   =   1/v∫_v_eτ  dv .

The temperature gradient at the border between the two parts of the boulder can be estimated as the temperature difference divided by the distance, (τ_w-τ_e)/l_ew, 
with l_ew being the typical distance between the eastern and western parts of the boulder,
i.e. roughly the distance between the centers of the two parts.
To estimate the last term in the right-hand side of Eqn. (<ref>), τ can be substituted by its mean value, i.e. τ_w and τ_e for the western and the eastern parts of the boulder respectively.[When interpreted literally, the assumptions of a constant temperature gradient in the body of the boulder and of a constant temperature on the two parts of its surface might seem to contradict each other. Still, these assumptions should provide an acceptable estimate for the corresponding terms in Eqn. (<ref>), and thus finally lead to a reasonable estimate for TYORP. 
It is possible to construct a more sophisticated model for the temperature distribution inside the boulder, but it will lead to a more complicated mathematics and a more obscure physics, while its accuracy will be anyway largely negated by the assumptions I am going to make below.]
Therefore, with these simplifications in place, Eqn. (<ref>) for the two parts of the boulder assumes the following form:

    θ^2 v τ_w'   =   s_st/l_ew(τ_e-τ_w)+s_spi_w-s_spτ_w^4 , 
    θ^2 v τ_e'   =   s_st/l_ew(τ_w-τ_e)+s_spi_e-s_spτ_e^4 .

Here i_w and i_e denote the dimensionless solar energy fluxes, defined as

    i_w   =   1/Φ s_sp∫_s_sp wI_ids_i , 
    
     i_e   =   1/Φ s_sp∫_s_sp eI_ids_i .


I decompose i_w(ϕ) and i_e(ϕ) into a Fourier series,
and disregard all of the terms except for the zeroth and first order ones.
This simplification can alter the final result, but still can serve as an estimate.
Even with a modified illumination function, it is still a valid physical problem,
whose solution must still bear the basic properties of the tangential YORP.
Thus, for insolation I substitute

    i_w=Ccosϕ+Ssinϕ+τ_0^4  , 
    
     i_e=Ccosϕ-Ssinϕ+τ_0^4  ,

where C, S, and τ_0 are constants.
Here, τ_0 has a physical meaning of the dimensionless temperature of the boulder, for which the emitted power equals the time-averaged absorbed power. Usually τ_0 is close to the mean temperature of the boulder.
When writing the same τ_0 in both equations, the same coefficients for cosine and opposite coefficients for sine,
I took into account that in the morning the eastern part of a symmetric boulder is illuminated in exactly the same manner
as the western part is in the evening, so that i_w(ϕ)=i_e(-ϕ).

Next, I introduce the coefficients

    a=v/s_spl ,    b=s_stl/s_spl_ew ,

where l is some typical boulder size.
Both coefficients a and b depend solely on the boulder shape, while the dependence on size only enters through l.

Finally, substituting Eqs. (<ref>) and (<ref>) into Eqn. (<ref>), I obtain

    θ^2 al τ_w'   =   b/l (τ_e-τ_w)+Ccosϕ+Ssinϕ+τ_0^4-τ_w^4 ,
    θ^2 al τ_e'   =   b/l (τ_w-τ_e)+Ccosϕ-Ssinϕ+τ_0^4-τ_e^4 .

Although this system is much simpler than the exact partial differential equation describing the heat conduction,
it still cannot be exactly solved analytically because of the nonlinearity τ^4.
Therefore, I aim to construct its approximate analytic solution, which I do in the following subsection.



 §.§ Approximate solution of the heat conduction equations

I am looking for the solution of Eqs. (<ref>) in the form of a series in terms of C and S,

    τ_w   =   τ_w0+τ_w1+τ_w2+... , 
    τ_e   =   τ_e0+τ_e1+τ_e2+... ,

where τ_w0 and τ_e0 are independent of C and S,
τ_w1 and τ_e1 are proportional to the first powers of C and S,
τ_w2 and τ_e2 are proportional to their second powers, and so on.
As I am looking for a relaxed periodic solution, all the terms have to be periodic in ϕ with a period of 2π.
I am going to account for the contribution to TYORP of only the first three terms, and to disregard higher-order terms.
This is perfectly justified if C,S≪τ_0^4, but usually C and S are only slightly less than τ_0 (compare with Table <ref>).
It implies that the decomposition in the form of Eqn. (<ref>) converges slowly, if at all.
Still, making C and S a factor of few smaller would make for a good convergence,
and all the further analysis would be justified.
Then I can assume that the extrapolation of the approximate formulas into the domain C∼ S∼τ_0^4 
must give a reasonable order-of-magnitude estimate for TYORP.

Note, that Eqn. (<ref>) is a Taylor series in terms of C and S, 
in contrast to Eqn. (<ref>), which is a Fourier series in terms ϕ.
I hold Fourier terms in Eqn. (<ref>) up to the first order and Taylor terms in Eqn. (<ref>) up to the second order, in order to construct the minimal model for TYORP.
With merely zeroth-order terms in either Eqn. (<ref>) or Eqn. (<ref>), TYORP would vanish.
Therefore, retaining first-order terms in both expansions is absolutely necessary. 
Retaining also the second-order Taylor term in Eqn. (<ref>) is motivated by the fact that its contribution is of the same order as the contribution of the first-order Taylor term (as we will see later from Eqs. (<ref>) and (<ref>)).
Taking more terms in either the Fourier decomposition Eqn. (<ref>) or the Taylor decomposition Eqn. (<ref>) can make the solution more precise, but at the cost of increasing complexity of the problem and of the final expression.
Moreover, treating even an infinite number of terms in both the decompositions will not make the solution exact, because the initial Equation (<ref>) is already an approximation.
This can turn any attempt to go beyond the minimalistic model in Eqs. (<ref>) and (<ref>) into an overkill.

Thus I substitute Eqn. (<ref>) into Eqn. (<ref>),
and equate the terms of the same order. 
In the highest (zeroth) order I get

    θ^2 al τ_w0'   =   b/l (τ_e0-τ_w0)+τ_0^4-τ_w0^4 , 
    θ^2 al τ_e0'   =   b/l (τ_w0-τ_e0)+τ_0^4-τ_e0^4 .

The periodic solution of this equation is τ_w0=τ_e0=τ_0.

I substitute the obtained τ_w0 and τ_e0 back into Eqs. (<ref>) 
and write the terms of the first order, which are linear in terms of C and S:

    θ^2 al τ_w1'   =   b/l (τ_e1-τ_w1)+Ccosϕ+Ssinϕ-4τ_0^3 τ_w1 ,
    θ^2 al τ_e1'   =   b/l (τ_w1-τ_e1)+Ccosϕ-Ssinϕ-4τ_0^3 τ_e1 .

This is a system of linear differential equations with a sinusoidal inhomogeneity,
whose periodic solution can be found in the form

    τ_w1   =   C_w1cosϕ+S_w1sinϕ ,
    τ_e1   =   C_e1cosϕ+S_e1sinϕ .

By substituting Eqn. (<ref>) into Eqn. (<ref>) and equating the coefficients in front of the sines and cosines, I get

    C_w1   =   4τ_0^3C/16τ_0^6+a^2l^2θ^4-al^3θ^2S/4(b+2lτ_0^3)^2+a^2l^4θ^4 ,
    
     S_w1   =   alθ^2C/16τ_0^6+a^2l^2θ^4+2l(b+2lτ_0^3)S/4(b+2lτ_0^3)^2+a^2l^4θ^4 ,
    
     C_e1   =   4τ_0^3C/16τ_0^6+a^2l^2θ^4+al^3θ^2S/4(b+2lτ_0^3)^2+a^2l^4θ^4 ,
    
     S_e1   =   alθ^2C/16τ_0^6+a^2l^2θ^4-2l(b+2lτ_0^3)S/4(b+2lτ_0^3)^2+a^2l^4θ^4 .


Finally, I write down the second order terms of Eqs. (<ref>),
which are quadratic in terms of C and S:

    θ^2 al τ_w2'   =   b/l (τ_e2-τ_w2)-6τ_0^2 τ_w1^2-4τ_0^3 τ_w2 , 
    θ^2 al τ_e2'   =   b/l (τ_w2-τ_e2)-6τ_0^2 τ_e1^2-4τ_0^3 τ_e2 .

I subtract Eqs. (<ref>) from each other, and average the result.
The left-hand side averages to 0, as τ_w2 and τ_e2 are periodic, and I am left with

    ⟨τ_w2⟩-⟨τ_e2⟩
     =-3lτ_0^2/b+2lτ_0^3 (⟨τ_w1^2⟩-⟨τ_e1^2⟩) .

One does not need to find the exact expressions for τ_w2 and τ_e2:
as will be seen in the next subsection, Eqn. (<ref>) suffices to compute TYORP
in the second order in terms of C and S.



 §.§ Computation of TYORP

A heated surface emits light and experiences the recoil pressure

    P=2/3cϵσ T^4=(1-A)Φ/c2/3τ^4 .

Here c is the speed of light, and the coefficient 2/3 corresponds to the light emission in accordance with Lambert's law.

To get the force experienced by the boulder in the eastward direction, I integrate this pressure over the boulder's surface, and take the F_1 component of the total force, assuming that the x_1 axis is directed from west to east:

    F_1=-∫_S_spP dS_1 .


I average this force over time and nondimensionalize it by dividing it by (1-A)Φ S_proj/c,
where S_proj is the horizontal projected area of the boulder.
Thus I get the dimensionless TYORP pressure,

    p=-2/31/s_proj∫_s_sp⟨τ^4⟩ ds_1 .

Now let us recall the assumption that the boulder is separated into two parts, 
with the dimensionless temperatures τ_w and τ_e,
and the two parts are symmetric, with the same surface area of the western and eastern parts, 
s_sp w=s_sp e.
Assuming ds_1 to always be negative at s_sp w and always positive at s_sp e, I get

    p   =   -2/31/s_proj(
    ⟨τ_w^4⟩∫_s_sp wds_1
    +⟨τ_e^4 ⟩∫_s_sp eds_1)
       =   2/31/s_proj(
    ⟨τ_w^4⟩∫_s_sp wd|s_1|
    -⟨τ_e^4 ⟩∫_s_sp ed|s_1|)
       =   2/3⟨τ_w^4-τ_e^4⟩1/s_proj∫_s_sp wd|s_1|
       =   2/3s_ew/s_proj⟨τ_w^4-τ_e^4⟩ .

Here s_ew=∫_s_sp wd|s_1| denotes the area of the boulder projected onto the vertical meridianal plane. It is convenient to introduce a new coefficient,

    n=2s_ew/s_proj .

This coefficient n characterizes steepness of the surface. If s_proj is understood as the horizontal projected area of the boulder, then n=⟨tanα_ew⟩ is the mean tangent of the slope of the surface in the east-west direction.
It is zero for a flat surface, unity for a surface with 45^∘ slopes, and bigger for even steeper slopes. 
Now, substituting Eqn. (<ref>) into Eqn. (<ref>) 
and using the decomposition of τ from Eqn. (<ref>), results into

    p   =   1/3n⟨τ_w^4-τ_e^4 ⟩
       =   1/3n⟨τ_0^4+4τ_0^3τ_w1+4τ_0^3τ_w2+6τ_0^2τ_w1^2+... 
          -τ_0^4-4τ_0^3τ_e1-4τ_0^3τ_e2-6τ_0^2τ_e1^2-... ⟩ .

As τ_w1 and τ_e1 are sinusoidal, their means are 0.
The means of τ_w2 and τ_e2 are obtained from Eqn. (<ref>). 
Thus Eqn. (<ref>) results in

    p=2nbτ_0^2/b+2lτ_0^3 (⟨τ_w1^2⟩-⟨τ_e1^2⟩) .

Using Eqs. (<ref>) for τ_w1 and τ_e1 leads to

    ⟨τ_w1^2⟩-⟨τ_e1^2⟩=1/2(C_w1^2+S_w1^2-C_e1^2-S_e1^2) .

Finally, substituting this expression into Eqn. (<ref>) and using Eqn. (<ref>), I obtain

    p=8nab^2τ_0^2l^2θ^2CS/(b+2lτ_0^3)(16τ_0^6+a^2l^2θ^4)(4(b+2lτ_0^3)^2+a^2l^4θ^4) .


This is the analytic expression of TYORP that I was looking for. It allows one to approximately estimate TYORP using the shape and thermal properties of the boulder.
In the following section, I apply it to several different boulder shapes.



§ APPLICATION OF ANALYTIC TYORP EXPRESSION TO BOULDERS OF DIFFERENT SHAPES



 §.§ One-dimensional wall

Let us start testing this analytic expression with the simplest and historically first model of TYORP:
one-dimensional heat conductivity in a wall by <cit.>.
The high long wall stands on regolith, going from the north to the south, and is illuminated by the sun alternatively from the east or the west (panel (a) in Figure <ref>).
Given that the dimensionless thickness of the wall is d and considering its patch of dimensionless surface area s, I separate the wall into the eastern and the western slabs with the thicknesses d/2. Then the volume of each slab is v=sd/2 and the typical distance between the slabs is l_ew=d/2. Thus I go on, filling in the geometric properties in the first row of Table <ref>, until I get a, b and n. (Note though, that for s_proj I take not a horizontal projected surface area, but a vertical one, which is the only meaningful definition for a very high vertical wall. It also renders n different from its conventional meaning of ⟨tanα_ew⟩.)

To compute the last three columns in Table <ref>, I need the insolation function, for which I take

    i_w(ϕ)=2sinϕ,  0<ϕ<π/2.

The coefficient 2 is due to the assumed mirror reflection of light by the regolith.
Decomposing i_w(ϕ) into the Fourier series, I finish filling in the first row of Table <ref>.









Now I substitute the obtained coefficients a, b, n, τ_0, C and S into Eqn. <ref>, and in Figures <ref> and <ref> compare the resulting analytic expression for TYORP pressure p with the numeric simulations by <cit.> and <cit.>.
The agreement is almost too good, given the approximations made while constructing the theory.
For the values of l and θ that correspond to big p, the value of p is predicted with an accuracy of a few tens per cent (Figure <ref>), although far from the maximal p the accuracy is worse.
The area where p is big, is very similar in the analytic theory and in the simulations (Figure <ref>).
The qualitative behavior of p as a function of l and θ is also reproduced correctly.
This is as good as one could expect from a simple estimate done by such an approximate theory.

A few warnings are still to be made.
In Figure <ref>, I disregard negative simulated values of p, observed by <cit.> for big l and θ.
<cit.> suggested that these negative values could be numeric artifacts,
and this was later confirmed by <cit.>, who only obtained positive p in all their simulations.
Furthermore, <cit.> proved that mirror reflection of light by the regolith is a bad approximation (worse, in fact, than disregarding reflected light completely).
Still, this approximation enters both the theory and the simulations, 
so they should be wrong in the same way, and thus agree with each other, which they do.



 §.§ Spherical boulder

The second model I want to study is a particular case from <cit.>:
a spherical boulder of radius r lying half buried in the regolith on the equator of an asteroid (panel (b) in Figure <ref>).
I separate the boulder into the eastern and the western hemispheres,
and list all the necessary geometric properties in Table <ref>.
For the typical distance between the hemispheres l_ew I take r, which is half the distance between their most remote points. Measuring the distance between their centers of mass instead would result into 3/4r, and an argument could be also made for using this value instead, although the difference should not matter much given the crudeness of all the previously made assumptions.
For s_proj I take the horizontal projection of the boulder, π r^2.
This is different from <cit.>, who considered a regular array of boulders and used for s_proj the surface area per boulder. Therefore, I re-normalize their results. Note that the choice of s_proj does not influence the observable physical values: s_proj enters the definition of p inversely, but also enters the transformation factor from p to F directly, and thus cancels out.

For the insolation I try two possibilities.
First, assuming mirror reflection of light by the regolith, I get (see Figure <ref> for explanation)

    i_w(ϕ)=1/2(1+sinϕ),  -π/2<ϕ<π/2.

Second, assuming full absorption of light by the regolith, and thus accounting for only the direct solar irradiation, I get

    i_w(ϕ)=1/2(1+sinϕ),  if   -π/2<ϕ<0, 
    
     i_w(ϕ)=1/2(cosϕ+sinϕ),  if   0<ϕ<π/2.
 
I treat the two cases separately in the second and third rows of Table <ref>.



In Figure <ref>, I compare the numeric simulations of <cit.> (solid lines)
with the two analytic models.
Qualitatively both of them work well, but the model including only the direct light has a better quantitative agreement with the simulations than the model with mirror reflection (similarly to findings of <cit.> for a wall).
In Figure <ref>, I compare the better analytical model with the numeric simulation.
One must also keep in mind that the simulations of <cit.> were conducted for an array of boulders partially shadowing each other.
Subtracting the shadowing should increase TYORP by a few tens per cent, as seen in the middle right panel of Figure 4 in <cit.>.
This could somewhat improve the agreement between the presented analytical model and the simulations.







 §.§ Wave in the regolith


Now, having tested the model for already studied cases, let us move to an as-yet unknown terrain.
Consider a regular array of hills and dales, with flat slopes angled at α, with a dimensionless wavelength λ,
all positioned on the equator in the north-south direction (panel (c) in Figure <ref>).
The model is intended as a proxy for small bumps and pits on the regolith, whose contribution to TYORP has never yet been considered.
What I am going to do here is but a crude estimate, yet to be tested by thorough numeric simulations.
Unlike for a boulder with a well-specified volume, now the volume in which heat conduction occurs has no well-defined borders, and estimating geometric properties for Table <ref> gets more complicated. Obviously, considering an infinitely deep volume below the surface is unsatisfactory, as an infinitely large volume will have infinite heat capacity, zero temperature oscillations, and thus zero TYORP.

Differently directed slopes cause uneven heating of the surface, and thus temperature variations in the east-west direction.
The principal Fourier harmonics of this horizontal temperature oscillation has the wavelength λ, and the corresponding wave vector k=2π/λ. The temperature oscillation will dampen with depth, and one can estimate the penetration depth of these temperature oscillations to be of the order of 1/k=λ/2π.
Thus I consider heat conduction in a slab of depth λ/2π and surface area s.
Although for λ or θ very different from unity the typical depth can be different and determined by other physical processes,
and thus the model can produce a bigger error than in the previous two cases,
even now it must at least give a reasonable order-of-magnitude estimate of TYORP for λ∼θ∼ 1.
The slab under consideration borders two other slabs, to the east and to the west.
These two slabs have the same temperature, and one can account for having two of them by just doubling their area of contact, s_st.

Having said all this, I fill in the last row in Table <ref>.
I assume that α is small, thus substituting cosα≈ 1, tanα≈sinα≈α.
For the illumination function I take

    i_w(ϕ)=cos(ϕ-α),  -π/2+α<ϕ<π/2-α.

When Fourier-decomposing i_w(ϕ), I also keep only the principal terms in terms of α.

Substituting the last row from Table <ref> into Eqn. <ref>, I get Figure <ref>.
The figure was produced assuming slope angles α=0.1. For other slope angles the TYORP pressure will scale as α^2. 
The maximal TYORP for regolith in Figure <ref> is about 30 times smaller than the maximal TYORP for spherical boulders in Figure <ref>.
Therefore, at first sight, it may seem that regolith is unimportant.
Still, one must bear in mind several reservations.
Firstly, non-cracked boulders have big thermal parameters, θ=10→ 100, while for regolith θ∼ 1 <cit.>.
This will substantially suppress TYORP for boulders, while leaving TYORP for regolith close to its maximum.
Secondly, a much bigger fraction of surface is covered by regolith than is by boulders.
Lastly, the typical slopes α for asteroids are utterly unknown, and altering α can drastically decrease or increase TYORP.
The length scales I am referring to are millimeters to centimeters, which is far below the resolution of even in situ observations of asteroids.

Summing up, one should acknowledge that regolith may have a substantial contribution to TYORP, which deserves a more detailed analysis than this rough estimate. Still, such an analysis would greatly divert us from the main focus of this paper.





§ INTEGRATION OF TYORP OVER ALL BOULDER SIZES



 §.§ Derivation of the integral for the overall TYORP


All the previous analysis was done for boulders of only one particular size.
Let us now integrate TYORP over all sizes of boulders on the asteroid,
thus generalizing the numeric results obtained by <cit.> for 25143 Itokawa.
For the size distribution of boulders, I take a power law: it is simple, reasonably precise, and the most widely used distribution. 
The interested reader is referred to <cit.> for an excellent review of the available measurements of the boulder size distribution for different asteroids.

Consider a power-law distribution of boulder sizes,

    dN/dl=N_0 l^-γ ,

where dN is the number of boulders with sizes between l and l+dl lying on the entire surface area of the asteroid s_ast, while N_0 and γ are constant.
For us it is more convenient to re-write this distribution not in terms of the number of boulders, but in terms of the area covered by them:

    dn/dl=n_0 l^2-γ .

Here dn=dS/S_ast is the fraction of the surface area covered by boulders with sizes between l and l+dl, while n_0=N_0 s_proj/s_astl^2 is a new constant. 
(Remember that for geometrically similar boulders of different sizes s_proj∝ l^2.)

Of course, a power law in the form of Eqs. (<ref>) or (<ref>) can only be an approximation.
This distribution would necessarily diverge for either small or big boulders, and would necessarily eventually start disagreeing with the data for both small and big boulders. Still, I assume that deviations from the power law happen only for such small and such big boulders
that their contribution to TYORP is negligible anyway.

Let us determine the total dimensionless TYORP drag p_tot experienced by the surface as the dimensional TYORP drag force divided by (1-A)Φ S/c.
Then p_tot can be obtained from p(l) by integrating it over fractions of the surface area dn covered by boulders of each size:

    p_tot=∫ p(l) dn=∫_0^∞ n_0 p(l)l^2-γ dl

Given asymptotics of Eqn. (<ref>) p∝ l^2 for l→ 0 and p∝ l^-5 for l→∞, this integral converges for γ∈(-2;5). Typical power indices for small boulders indeed lie in this range <cit.>. For γ outside this interval, finiteness of the integral has to be provided by setting the minimal or the maximal size of boulders, starting from which the power-law size distribution breaks.

From Eqn. (<ref>) one can see the problem: the rational expression for p provided by Eqn. (<ref>) does not look frightening only up until the moment one tries to integrate it over l.
Having spent many pleasant hours doing contour integration in the complex plane, I must acknowledge that the final result is too complicated to be of any practical use despite it being a closed form algebraic expression.
Moreover, when written in terms of real variables, it is far too lengthy to be accommodated into this subsection.

To get a more practical result from Eqn. (<ref>), p should be transformed into a simpler form.
Still, one can not neglect either term in Eqn. (<ref>): a, b, τ_0, l and θ are all of the order of unity, and thus all comparable with each other.
Therefore, I choose not to simplify the equation, but to brute-forcedly interpolate it.
This will cause some loss of accuracy, but as we will see, this loss is not much bigger than the errors already caused by the previous simplifications.



 §.§ Approximate expression for TYORP


In this subsection I do a detour to find the best fit for p as a function of l and θ as given by Eqn. (<ref>).
In Figures <ref> and <ref> one can see the three-dimensional surface p(l,θ) from above, and notice that in the l-θ plane lines of constant p are roughly elliptical and roughly similar to each other.
In Figures <ref> and <ref> one can see that vertical cross-sections of this three-dimensional surface look similar to Gaussian functions of the same width.
These observations allow us to guess that a decent approximation to p(l,θ) will be given by a two-dimensional log-normal function,

    p=p_0e^A_llnl+A_θlnθ+1/2A_llln^2l+A_lθlnllnθ+1/2A_θθln^2θ.


There are several ways to choose the coefficients in this fitting function, and in Figure <ref> I explore different ways of fitting.
The least squares fit to the numeric solution is presented with the orange lines,
the least squares fit to the approximate solution Eqn. <ref> is plotted in blue.
Both fits are constructed for the same range of r and θ as presented in Figure <ref>.
The purple lines are constructed by requiring Eqs. (<ref>) and (<ref>) to have the same values in the point l=θ=1, as well as the same first and second partial derivatives. For brevity, I call this approach the Taylor fit, as I have indeed requested equality of the zeroth, first and second order terms in the Taylor series for p(l,θ) and its fit.
The analytically computed coefficients of the Taylor fit are presented in Table <ref>.
Although cumbersome, they can still be instrumental to estimations of TYORP for complex-shaped bodies.

The exact numeric solution is overplotted with a red line, and the approximate solution given by Eqn. <ref> with a green line.
One can see that the least squares fit to the numeric result agrees with them very well, as do the two other fits agree with the approximate solution Eqn. <ref>. The Taylor fit to Eqn. <ref> is by construction perfect in the vicinity of the point l=θ=1, but far away from this point it also works quite well.

Coefficients of the three considered fittings are given in Table <ref>.

As a minor side result, in the same table I present α_0=tan(1/2arctan2A_lθ/A_ll-A_θθ), which determines orientation of the ellipsis with constant p for the log-normal distribution. In the l-θ plane decrease of p is slowest along the line θ∝ l^α_0. We see, that α_0 obtained from different models are consistent with each other, and roughly consistent with the elongation of the red and yellow areas in Figure <ref>.

Now, having obtained the coefficients of the log-normal fit, I am ready to finish computing the integral for the overall TYORP.











 §.§ Computing the integral for the overall TYORP



Returning to the main focus of this section,
I still need to integrate the simplified expression for TYORP over all boulder sizes.
Substituting Eqn. (<ref>) into Eqn. (<ref>) and computing the Gaussian integral, I get

    p_tot=n_0μe^-(lnθ-lnθ_0)^2/ν^2,

where three new constants have been introduced:

    μ   =   p_0√(2π/-A_ll)e^2A_θ A_lθ(A_l-γ+3)-A_θθ(A_l-γ+3)^2-A_θ^2A_ll/2(A_θθA_ll-A_lθ^2), 
    ν   =   √(-2A_ll/A_θθA_ll-A_lθ^2), 
    lnθ_0   =   A_lθ(A_l-γ+3)-A_θ A_ll/A_θθA_ll-A_lθ^2.


The mathematical expressions are again somewhat lengthy, but they are worth the agreement between the different estimates of p_tot we see in Figure <ref>.
Even the substantial discrepancies between the numeric and the analytic solutions observed in Figure <ref> are mostly eliminated by integration in Figure <ref>.
The resulting discrepancies between the exact results and the different approximations lie within about 20% of the maximal value on the plot,
thus certifying that both of our approximations are good for estimating the TYORP effect.
The constants Eqn. (<ref>) are also listed in Table <ref>, and one can see that they are close to each other.

I now apply these results to asteroid 25 143 Itokawa, both to estimate the TYORP and to appease myself by certifying the uncertainties in physical properties of even such a well-studied asteroid cause a much bigger error in the TYORP acceleration than the simplifications underlying my analytic model.
The assumed physical properties of boulders taken from <cit.> are listed in Table <ref>.
Following <cit.>, I use two different values of the heat conductivity of rock: bigger κ_1 presumably corresponds to solid rock, while smaller κ_2 could correspond to cracked rock, which is assumed to be widely present on atmosphereless bodies <cit.>.
Respectively, I get two different thermal parameters θ_1=13 and θ_2=4.1.
I use the size distribution of boulders from <cit.>, dN=(14± 9)· 10^3 L^-3.0± 0.2dL, where dL is in meters.
For simplicity, I choose the power index γ=3. In this case n_0 is independent of L_cond. 
Nondimensionalizing L, transforming from the number of boulders dN into their relative surface area dn, assuming roughly circular projections of boulders, and taking for the surface area of Itokawa 3.93· 10^5 m^2 <cit.>, I get n_0=0.028± 0.018.
From Figure <ref>, I find corresponding pressures P(θ_1)≈ 0.001 and P(θ_2)≈ 0.005.
Estimating the dimensionless TYORP of the whole asteroid as τ_z≈ 9sT(θ) <cit.>, I get for solid rocks τ_z1=0.00025± 0.00015 and for cracked rocks τ_z2=0.0013± 0.0007.
This value should be compared to the difference Δτ_z= 0.003± 0.002 <cit.> between the observed accelerating torque acting on Itokawa and the predicted normal YORP torque, as one may expect Δτ_z to be equal to Itokawa's TYORP. One can see that the cracked rock gives a value consistent with Δτ_z, while the solid rock gives a value an order of magnitude smaller than Δτ_z.
As a side remark, I note that the tens-of-per-cent uncertainty of the theoretical model is much smaller than an order of magnitude uncertainty caused by the uncertainty of the physical properties of the asteroid.
The most important takeaway message from this estimate for 25 143 Itokawa is that for realistic properties of boulders and their realistic shape distribution, Eqn. (<ref>) predicts a significant effect, comparable to the normal YORP.





§ DISCUSSION


Reviewing the derivation of the analytic expression of TYORP (Eqn. (<ref>)), one can single out three major ingredients sufficient for TYORP to appear:

1. Heating pattern with a day–night asymmetry (C≠ 0) and an east–west asymmetry (S≠ 0).

2. Thermal inertia (0<θ<∞, 0<a<∞) and thermal conduction (0<b<∞).

3. Non-linearity of the heat emission law.

The non-linearity enters the equations for TYORP in two related ways: through the boundary condition of the heat conduction equation (Eqn. (<ref>)) and through the definition of TYORP pressure (Eqn. (<ref>)).
The former's contribution is negative, the latter's contribution is positive, and the latter wins, as can be seen from the substitution of Eqn. (<ref>) into Eqn. (<ref>).

Acknowledging the generality of conditions 1–3, one must conclude that TYORP is a very general feature, which is not bound to any particular boulder shape or illumination pattern. More or less any boulder shapes can be substituted into Eqn. (<ref>), while Eqn. (<ref>) presents a most generic illumination function, which is still capable of causing TYORP of a realistic magnitude, although boulders with steeper slopes (bigger n) and illuminated in a more asymmetric way (bigger S) are subject to systematically bigger TYORP, as can be seen from Eqn. (<ref>).

In my model, C is always positive, as the boulder is illuminated in the day, not at night.
S is usually positive, as all realistic boulders are more illuminated from the east before noon,
and more illuminated from the west after noon.
Thus all terms in Eqn. (<ref>) are positive, and the resulting TYORP derived from the analytic model should accelerate the rotation of the asteroid under any realistic conditions, not decelerate it.
This fact has already been observed for different geometries of boulders by <cit.> and <cit.> in their numerical simulations.

This holds not only for convex structures like boulders or mounds, but also for concave structures, like grooves or pits on the surface of the asteroid.
In the latter case the mechanism is the same as usual: east-facing slopes are better illuminated in the morning, west-facing slopes in the afternoon, and the heat conduction occurs between them. 
The problem is again approximately described by Eqs. (<ref>) and (<ref>) with all constants being positive, and the solution again leads to Eqn. (<ref>).
It means that the west-facing slopes again emit on average more light, push the surface more in the eastern direction, and again accelerate the rotation of the asteroid.

Thus most imaginable surface structures should be to a certain extent subject to positive TYORP: crater rims, crater pits, grooves and cracks on the surface, small mounds and summits, or just wavy patterns of regolith, rough surfaces of stones, lone boulders and boulder fields... A heap of pebbles can contribute to TYORP either as a single body, or via individual pebbles, or via individual pocks and wedges of individual pebbles, and which of these descriptions is the most appropriate depends on the rotation rate of the asteroid, the heat conductivity of rock, as well as other parameters. Singling out the most important contributors to TYORP and accounting for them all is an important task for future research.



§ CONCLUSIONS


I have constructed an approximate fully analytic model of TYORP, in which I used an averaging procedure to transform a partial differential equation into an ordinary differential equation, then used perturbation theory to transform the latter into algebraic equations, and finally arrived at an algebraic expression for TYORP.
Along the way the shape of the boulder was boiled down to just a few geometric coefficients, and the illumination pattern – to its few Fourier harmonics.
Despite the simplicity, the analytic expression for TYORP was found to be in good agreement with numeric simulations for a stone wall <cit.> and a spherical boulder <cit.>.

I used the analytic model to estimate the TYORP produced by the non-smoothness of regolith on an asteroid.
On the one hand, regolith should be much flatter than boulders, which diminishes its TYORP.
On the other hand, it has a more favorable thermal parameter θ and presumably covers a bigger fraction of the surface, which increases its TYORP contribution.
Given the utter uncertainty of the unevenness of regolith on different scales, it is hard to say whether or not its contribution is significant, but neglecting it now seems unsafe.

I then have integrated TYORP over boulders of all different sizes, assuming a power-law size distribution for boulders, and doing a two-dimensional log-normal fit to the analytic TYORP drag.
The additional error introduced by this fit was generally not bigger than the error already caused by the assumptions underlying the analytical model.
Moreover, the error substantially diminished after integration over the boulder sizes.
The resulting TYORP drag of the entire assembly of boulders appeared log-normal in terms of the thermal parameter θ, and concurrently in terms of rotation rate of the asteroid.
It means that even for a broad distribution of boulders over sizes, TYORP still has a relatively narrow peak, a corresponding rotation rate with the biggest TYORP, and a relatively fast decrease of TYORP for faster or slower rotation.
The integrated TYORP drags for the log-normal approximation, for the full analytic solution and for the numeric simulation were found to be in good agreement with each other.

The proposed analytic model demonstrates that TYORP should appear for most realistic shapes of boulders or regolith, and should	 be positive in most realistic cases.



§ ACKNOWLEDGEMENTS

I am very grateful to Daniel J. Scheeres for discussing with me the astronomical content of the problem, to Uliana Pyrohova for counseling me with mathematical issues, to Vlad Unukovych for reviewing the physics of the paper, and to Alexander Kostenko for improving the language style.
Without the help of either of them the paper would have been stuck at a certain point.

99
[Delbo et al.2014]delbo14
Delbo M., Libourel G., Wilkerson J., et al., 2014, Nature 508, 233 
[Demura et al.2006]demura06
Demura H., Kobayashi S., Nemoto E., et al., 2006, Science 312, 1347
[Golubov & Krugly2011]golubov11
Golubov O., Krugly Yu. N., 2011, Planetary Defence Conference, From Threat to Action, held 09-12 May 2011, in Bucharest, Romania 
[Golubov & Krugly2012]golubov12
Golubov O., Krugly Yu. N., Tangential component of the YORP effect, 2012, ApJL 752, L11
[Golubov et al.2014]golubov14
Golubov O., Scheeres D. J., Krugly Yu. N., A three-dimensional model of tangential YORP, 2014, ApJ 794, 22
[Murdoch et al.2015]asteroids4
Murdoch N., Sánchez P., Schwartz S. R., Miyamoto H., 2015, Asteroids IV, 767
[Ševeček et al.2015]sevecek15
Ševeček P., Brož M., Čapek D., & Ďurech, J., The thermal emission from boulders on (25143) Itokawa and general implications for the YORP effect, 2015, MNRAS 450, 2104
[Ševeček et al.2016]sevecek16
Ševeček, Golubov O., Scheeres D.J., Krugly Yu. N., Obliquity dependence of the tangential YORP, 2016, A&A 592, 115




