%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% A general linear theory of the Yarkovsky effects on spheres
% Revised by PF Dec.  30, 1998
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\documentstyle{l-aa}
\newcommand{\etal}{{\it et~al.}}
\newcommand{\etall}{{\it et~al. }}
%\renewcommand{\baselinestretch}{2}
%\onecolumn
\thesaurus{05.03.1; 07.13.2}
\title{A complete linear model for the Yarkovsky thermal force
 on spherical asteroid fragments}
\author{David Vokrouhlick\'y}
\institute{Institute of Astronomy, Charles University,
    V Hole\v{s}ovi\v{c}k\'ach 2, CZ-180 00 Prague 8, Czech Republic, and\\
           Observatoire de la C\^ote d'Azur, dept. CERGA, Av N. Copernic,
    F-06130 Grasse, France\\
    (Electronic address: vokrouhl@mbox.cesnet.cz)}
\date{\today}
\offprints{Prague address.}
\begin{document}
\maketitle
\markboth{D. Vokrouhlick\'y: Yarkovsky effect on spherical fragments...}{}
\begin{abstract}
A linear theory for heat conduction in a spherical, solid and
rotating body illuminated by solar radiation is developed.  The recoil
force due to the thermally re-emitted radiation by the surface of the
body is computed, including all the terms depending both on the body's
rotation frequency and the mean motion of its revolution about the
Sun. The present solution thus overcomes a drawback of the previous
approaches, which have been tailored separately either to the diurnal
or to the seasonal variant of the so--called Yarkovsky effect,
corresponding to different limiting cases of the current theory.  We
pay a special attention to compute the secular effects on the
semimajor axis of the body's orbit about the Sun. The results from the
general model coincide with those of the previous approaches to a high
level of accuracy, as the relative size of the additional ``mixed''
terms is smaller than $10^{-3}$ for plausible parameter choices.  This
confirms that the use of the simplified formul{\ae} is warranted in
the relevant Solar System applications.
\keywords
{Celestial mechanics -- Minor planets, asteroids -- Meteoroids}
\end{abstract}
 
\section{Introduction}
The so--called Yarkovsky effect, a recoil force due to thermal
radiation from anisotropically heated orbiting bodies, has recently
attracted a considerable attention in the frame of the studies on the
delivery of meteorites and the dynamics of small bodies in the Solar
System.  Specific issues for which the Yarkovsky effect is probably
relevant are: the cosmic--ray exposure ages of stony and iron
meteorites, which are much longer than the dynamical lifetimes of
particles delivered from the asteroid belt (Farinella et~al. 1998;
Hartmann et~al. 1998; Morbidelli \& Gladman 1998); the overabundance
of decameter--sized near--Earth objects (Rubincam 1995, 1998;
Vokrouhlick\'y \& Farinella 1998a); the dynamical evolution of
large ($>5$~km) main--belt asteroid fragments and their delivery to
Mars-- and Earth--crossing orbits (Farinella \& Vokrouhlick\'y 1998).
In all these cases, the Yarkovsky effect plays the role of a
dissipative mechanism, resulting into a significant long--term
mobility of the orbital semimajor axis and a complex interaction with
resonances.
 
To assess the relevance of the Yarkovsky effect in Solar System
dynamics one needs, as a first step, to develop a reliable physical
model of the thermal processes occurring within solid, spinning and
orbiting bodies.  A significant amount of work has been performed on
this problem in recent years, after Rubincam (1995) resurrected the
interest in the dynamical consequences of these thermal effects.  Most
importantly, Rubincam (1995, 1998) and Farinella et~al.  (1998)
recognized the existence of two distinct variants of the Yarkovsky
effect: a ``diurnal'' variant depending on the rotation frequency
of the body around its instantaneous spin axis ($\omega_{\rm rot}$),
and a ``seasonal'' variant depending on the mean motion frequency
of the body around the Sun ($\omega_{\rm rev}$).
 
Technically speaking, the diurnal variant is obtained when one
entirely neglects the orbital motion around the Sun (see e.g.
Vokrouhlick\'y 1998a,b), whereas in dealing with the seasonal variant
one {\it a priori} averages all relevant quantities over the
(assumedly) fast rotation of the body (e.g.  Rubincam 1995, 1998;
Vokrouhlick\'y \& Farinella 1998b).  This classification is meaningful
and useful, since the two variants of the Yarkovsky effect result in
qualitatively different long--term changes of the semimajor axis.  The
diurnal version is maximum at zero obliquity and can lead either to
semimajor axis decrease or increase, depending on the sense of
rotation; on the contrary, the seasonal version is maximum at
$90^\circ$ obliquity and can only result in orbital decay (e.g.
Rubincam 1995, 1998; Farinella et~al. 1998; Hartmann et~al. 1998).  At
the essence, however, the two variants of the Yarkovsky effect are
just two different limiting cases of a single physical mechanism,
i.e., the recoil force associated to thermal radiation from a body
having an anisotropic temperature distribution on its surface.  As
their names imply, the diurnal and seasonal variants correspond to
different periodicities and geometries of the external illumination on
the body's surface.  From this perspective, it seems desirable to
develop a unified, self--consistent theory for the Yarkovsky effect,
including at the same time both the diurnal and the seasonal
periodicities, such that the two classical variants can be derived
computing suitable mathematical limits.
 
Although the classical variants of the Yarkovsky effect are present as
particular limiting cases, the unified theory inevitably will contain
additional, ``mixed'' terms, depending on both the relevant
frequencies $\omega_{\rm rot}$ and $\omega_{\rm rev}$.  This
conclusion holds even in the frame of a linear theory for the
temperature changes, such as that developed in the following Sections.
Thus, the major novelty of this paper consists of the derivation of
these ``mixed'' (or ``diurnal--seasonal'') terms.  Specifically, we
shall show that the ``diurnal'' variant does not exist as an effect
depending on the rotation frequency alone, but inevitably contains a
linear combination of the two frequencies.  As expected, in the limit
of a very rapid spin rate this doublet tends to merge into a single
line, depending just on $\omega_{\rm rot}$.
 
Then, we shall assess the contribution of the new terms to the secular
changes in the semimajor axis of the body's orbital motion.  As noted
above, such changes probably play an important role in several
problems of astronomical interest, and the quantitative results
obtained so far have always been computed as a simple superposition of
the diurnal and seasonal effects (e.g.  Farinella \& Vokrouhlick\'y
1998), neglecting any possible ``mixed'' effects.
 
To make the calculations as simple as possible, we shall make three
simplifying assumptions: (i) a circular orbit around the Sun; (ii) a
spherical shape of the body; and (iii) a commensurability between the
rotation and revolution periods.  In particular, we shall introduce a
parameter $m = \omega_{\rm rot}/\omega_{\rm rev}$, and we shall assume
that $m$ is an integer number.  However, we stress that while the
first two assumptions correspond to physical simplifications, the third
one is just a suitable mathematical step to simplify the derivation of
our results, and that this assumption can be easily removed by the
technique used by Farinella \& Vokrouhlick\'y (1996).  Therefore, our
final results will be valid for any (real) value of the parameter $m$.
 
\section{Theory}
\subsection{Formulation of the problem}
 
Hereafter, we use the mathematical approach and the notations
introduced by Vokrouhlick\'y (1998a,b).  We refer to those papers
for a historical background and a more detailed discussion of
thermal physics, while here we shall just provide a few general
concepts and definitions required to derive a unified solution for the
thermal Yarkovsky force on a spherical body.
 
Since we plan to remain in the framework of a linear theory to
describe the thermal response of an orbiting body to external heating,
we suppose that the temperature $T$ throughout the body is close to a
constant mean value $T = T_{\rm av}$, and therefore write $T = T_{\rm
av} + \Delta T$ ($\Delta T \ll T_{\rm av}$).  A suitable scaling of
the temperature $T$, as well as other variables, simplifies the
mathematical formulation of the problem.  The temperature $T$ will be
normalized by an auxiliary value $T_\star$ defined by: $\epsilon
\sigma T_\star^4 = \alpha {\cal E}_\star$.  Here, $\epsilon$ is the
thermal emissivity coefficient of the body, $\sigma$ the
Stefan--Boltzmann constant, $\alpha$ the optical absorption
coefficient and ${\cal E}_\star$ the solar radiation flux at the mean
distance along the orbit.  Since we shall restrict our analysis to the
case of a circular orbit, the scaled mean temperature reads $T'_{\rm
av} = T_{\rm av}/T_\star = 1/\sqrt{2}$ [see Vokrouhlick\'y 1998a; the
normalized quantities will be always denoted by a prime].  Similarly,
the radial coordinate $r$ measured from the center of the body to its
surface (at $r=R$) is to be scaled by the penetration depth of the
seasonal thermal wave $l_s = \sqrt{K/ \rho C \omega_{\rm rev}}$: $r'=
r/l_s$.  Here, $K$ is the thermal conductivity, $C$ the thermal
capacity, and $\rho$ the density of the material.  Finally, the time
$t$ will be represented by a complex quantity $\zeta = \exp(i
\lambda)$ [$\lambda=\omega_{\rm rev}(t-t_0)$], with $t_0$ being an
arbitrary time origin to be specified below (here $i = \sqrt{-1}$).
 
The temperature variation $\Delta T' = \Delta T/T_\star$
satisfies
the heat diffusion equation (see, e.g., Vokrouhlick\'y 1998a)
\begin{eqnarray}
 i \zeta {\partial \over \partial \zeta}\,& &\!\!\!\!\!\!\!\!\!
  \Delta T'(r';\theta,\phi;\zeta) = \label{one} \\&=&
   {1\over r'^2}\left\{{\partial \over \partial r'}
  \left(r'^2{\partial \over \partial r'}\right) + \Lambda\left(\theta,
  \phi\right)\right\}\, \Delta T'(r';\theta,\phi;\zeta) \; , \nonumber
\end{eqnarray}
with the operator $\Lambda(\theta,\phi)$ given by
\begin{equation}
 \Lambda\left(\theta,\phi\right) = {1\over \sin\theta}\left[{\partial
  \over \partial \theta}\left(\sin\theta{\partial \over \partial \theta}
  \right)+ {1\over \sin\theta}{\partial^2 \over \partial \phi^2}\right]
 \label{two}
\end{equation}
(here $\theta$ and $\phi$ are the usual spherical coordinates with the
pole defined by the body's rotation axis).
Eq. (\ref{one}) is supplemented by the linearized boundary
condition
\begin{equation}
 \sqrt{2} \Delta T'+ \Theta\left({\partial \Delta T'\over \partial r'}
  \right)_{R'} = \Delta {\cal E}'\; ,\label{three}
\end{equation}
with the seasonal thermal parameter $\Theta$ defined by
\begin{equation}
 \Theta = {\Gamma \sqrt{\omega_{\rm rev}}\over \epsilon\sigma
  T_\star^3} \label{four}
\end{equation}
(here $\Gamma = \sqrt{\rho C K}$ is the thermal inertia).  The
right--hand side term in Eq. (\ref{three}) is defined by ${\cal E}' =
\frac{1}{4} + \Delta {\cal E}'$, ${\cal E}'_{\rm av}=1/4$ being the
averaged irradiation of the fragment's surface.  The temperature
variation $\Delta T'$ is a function of the coordinates
$(r',\theta,\phi)$ inside the body and of time $\zeta$.
 
Like in Vokrouhlick\'y (1998a), we shall solve for $\Delta T'$ in a
rotating, body--fixed reference frame, with the $Z$-axis coinciding
with its spin vector.  At the reference time $t_0$, the $X$-axis of
this system points toward the radiation source (the Sun).  After
expressing the thermal force vector ${\bf f}$ in this reference frame
-- see Eqs.~(\ref{fifteen}) and (\ref{sixteen}) below -- we will
compute its projections in the orbit--related reference system and
obtain the force components appearing in the Gauss perturbation
equations.  The latter operation is simple vectorial algebra.
 
\subsection{The insolation term}
As discussed in detail by Vokrouhlick\'y (1998a), a particular
attention has to be paid to a suitable development of the irradiation
term $\Delta {\cal E}'$ on the right--hand side of Eq. (\ref{three}).
In general, we can use a series of spherical functions
\begin{eqnarray}
 \Delta {\cal E}' =\sum_{n\geq 1}\,\sum_{k=-n}^n a_{nk}(\zeta)\,
   Y_{nk}(\theta,\phi)\; , \label{five}
\end{eqnarray}
{} where only the dipole ($n=1$) part
\begin{eqnarray}
 a_{10}(\zeta) &=& \sqrt{{\pi \over 3}}\, \cos\theta_0(\zeta)\; ,
  \label{six} \\
 a_{1\pm 1}(\zeta) &=& \mp\sqrt{{\pi \over 6}}\, \sin\theta_0(\zeta)
 \, e^{\mp i \phi_0} \label{seven}
\end{eqnarray}
will be relevant. Here $\theta_0$ and $\phi_0$ are the solar
colatitude (measured
from the spin axis) and longitude in the body--fixed reference frame
(a
mathematical technique allowing one to obtain most elegantly these
results has been discussed in Vokrouhlick\'y 1998a). To
simplify the following algebra, the origin of the mean longitude
$\lambda$ (i.e. $t_0$) has been chosen so that
\begin{equation}
 \cos\theta_0 = -\sin\gamma \sin\lambda = {i\over 2}\sin\gamma
  \left(\zeta-\zeta^{-1}\right)\; , \label{eight}
\end{equation}
where $\gamma$ is the obliquity of the spin axis.
Eq.~(\ref{eight}) allows us to make immediately an
important observation, namely that the Fourier spectrum of the
$a_{10}(\zeta)$ coefficient contains only the revolution
frequency.
 
On the other hand, the Fourier development for the tesseral
coefficients $a_{1\pm 1}(\zeta)$ is more complicated. It contains
a linear combination of the rotation and revolution frequencies,
since
\begin{equation}
 \sin\theta_0\, e^{\pm i\phi_0} = \sin^2{\gamma\over 2}\; \zeta^{\mp (m+1)}
  + \cos^2 {\gamma\over 2}\; \zeta^{\mp (m-1)}\; . \label{nine}
\end{equation}
As an alternative way of obtaining these results, one may use
the insolation function development of Rubincam [1994; Eq.~(6)].
 
\subsection{Linear solution and thermal force evaluation}
The linearity of the system (\ref{one}) -- (\ref{three}) and the
development (\ref{five}) of the insolation term allows us to
make a
suitable decomposition of the temperature $\Delta T'$ into
multipole components
\begin{equation}
 \Delta T'(r';\theta,\phi;\zeta) = \sum_{n\geq 1}\,\sum_{k=-n}^n
  t'_{nk}(r';\zeta)\,Y_{nk}(\theta,\phi)  \; . \label{eleven}
\end{equation}
The coefficients $t'_{nk}(r';\zeta)$, weighing the different multipole
terms, satisfy the following system of decoupled differential
equations
\begin{eqnarray}
 i \zeta {\partial \over \partial \zeta}\, & & \!\!\!\!\!\!\!
  t'_{nk}(r';\zeta) =\label{twelve}\\
 &=&  {1\over r'^2}\left\{{\partial \over \partial r'}
  \left(r'^2{\partial \over \partial r'}\right) - n\left(n+1
  \right)\right\}\, t'_{nk}(r';\zeta)\; ,\nonumber
\end{eqnarray}
with the boundary constraints
\begin{equation}
 \sqrt{2}\, t'_{nk}(R';\zeta) + \Theta\left({\partial t'_{nk}\over
 \partial r'}\right)_{(R';\zeta)} = a_{nk}(\zeta)
 \label{thirteen}
\end{equation}
at the surface of the body. The regularity of $t'_{nk}$ at the central
position $r'=0$ is also assumed. In the next sections we shall obtain
a general solution of Eqs. (\ref{twelve}) and (\ref{thirteen}) for the
dipole ($n=1$) part of development (\ref{five}).
 
When determining the recoil force due to thermal
radiation, we assume -- in agreement with the boundary condition
(\ref{three}) -- that the isotropic Lambert's law holds as far as the
directional characteristics of the emission are concerned. Linearizing
the fourth power of the surface temperature as before, we obtain
\begin{equation}
 {\bf f}(\zeta) = -{2\sqrt{2}\over 3\pi} \alpha\Phi \int d\Omega\,
  \Delta T'(R';\theta,\phi;\zeta)\, {\bf n}\; ,\label{fourteen}
\end{equation}
which gives the thermal recoil force per unit of mass of the body.
Here,
$\Phi = ({\cal E}_\star \pi R^2/mc)$ is the usual radiation force
factor: $m$ is the body's mass, $c$ the velocity of light and ${\bf
n}$ the unit vector normal to the surface.
Given the multipole development (\ref{eleven}) of the
temperature variation $\Delta T'$, we easily obtain the
following formul{\ae} for the thermal force components
\begin{eqnarray}
 f_X(\zeta) + if_Y(\zeta) &=& -{8\over 3\sqrt{3\pi}}\, \alpha\Phi\,
  t'_{1-1}(R';\zeta)\; , \label{fifteen} \\
 f_Z(\zeta) &=& -{4\over 3}\sqrt{2\over 3\pi}\,\alpha\Phi\,
  t'_{10}(R';\zeta)\; , \label{sixteen}
\end{eqnarray}
{}confirming that only the dipole part ($n=1$) of the temperature
development (\ref{eleven}) is relevant at the level of a linear
theory. The force components
(\ref{fifteen}) and (\ref{sixteen}) are given in the co--rotating
body--fixed reference frame defined in Sec.~2.1.
 
\subsection{Solution for the seasonal component}
Hereafter, we give a solution for the ``zonal'' part $t'_{10}
(r';\zeta)$ of the temperature distribution, which yields the
along--spin component of the thermal force [see Eq. (\ref{sixteen})].
Given the simplicity
of the corresponding $a_{10}(\zeta)$ coefficient from
Eq. (\ref{eight}),
we may assume: $t'_{10}(r';\zeta) =\kappa_+(r')\zeta +\kappa_-(r')
\zeta^{-1}$. The radial functions $\kappa_\pm$ satisfy the equation
\begin{equation}
 \left[{d \over d r'}\left(r'^2{d \over d r'}\right) -\left(2
  \pm i r'^2\right)\right]\,\kappa_\pm(r') = 0 \; .
  \label{seventeen}
\end{equation}
whose solution is given by
\begin{equation}
 \kappa_\pm (r') = c^\pm  j_1(\sqrt{\mp i}r')\; .\label{eighteen}
\end{equation}
Here, $j_1(z)$ is the spherical Bessel function of order $1$.
Determining the proportionality factors $c^\pm$ in
Eq. (\ref{eighteen}) by the surface boundary
constraint (\ref{thirteen}), we obtain the following
expression for the zonal dipole coefficient $t'_{10}(r';\zeta)$
at the surface of the body ($r'=R'$):
\begin{equation}
 t'_{10}(R';\zeta) = -\sqrt{\pi\over 6}\, \sin\gamma \,
  {E_{R'} \sin(\lambda + \delta_{R'})\over 1+\chi}\; .
  \label{nineteen}
\end{equation}
Following the notations of Vokrouhlick\'y (1998a), we introduce
the amplitude $E_{R'}$ and the phase $\delta_{R'}$ by
\begin{equation}
 E_{R'}\exp\left(i\delta_{R'}\right) =
 {A(x)+iB(x)\over C(x)+iD(x)} \; , \label{twenty}
\end{equation}
with $x=\sqrt{2}R'$. The auxiliary functions
$A(x)$, $B(x)$, $C(x)$, $D(x)$ read
\begin{eqnarray}
 A(x) &=& -\left(x+2\right)- e^x\left[\left(x-2\right)\cos x - x
  \sin x\right]\; , \label{twentyone} \\
 B(x) &=& -x - e^x\left[x\cos x +\left(x-2\right)\sin x\right]\; ,
  \label{twentytwo} \\
 C(x) &=& A(x)+ {\chi\over 1+\chi}\times\label{twentythree} \\ & &
   \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \left\{3\left(x+2\right)+ e^x\left[
   3\left(x-2\right)\cos x+ x\left(x-3\right)\sin x\right]\right\}\;
   ,\nonumber \\
 D(x) &=& B(x)+ {\chi\over 1+\chi}\times\label{twentyfour}\\& &
  \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \left\{x\left(x+3\right) -e^x
  \left[x\left(x-3\right)\cos x- 3\left(x-2\right)\sin x\right]\right\}
  \; ,\nonumber
\end{eqnarray}
with the parameter $\chi$ in Eqs. (\ref{nineteen}),
(\ref{twentythree}) and (\ref{twentyfour}) defined
by: $\chi = \Theta/(\sqrt{2}R')$. Note that the same quantity had
been called $\lambda$ by Vokrouhlick\'y (1998a).
 
As we are mainly interested in the mean rate of change $(da/dt)$ of
the orbital semimajor axis $a$ due to Yarkovsky
effects, we look for the perturbations caused by the
along--spin thermal force component $f_Z$ in the along--track
direction -- $T_1$. Assuming a quasi--circular orbit we have
\begin{equation}
 T_1 = f_Z \sin\gamma \cos\lambda + {\cal O}(e)\; .
 \label{twentyfive}
\end{equation}
Averaging over one revolution of the body around the Sun (i.e.
one cycle of the mean longitude $\lambda$) we finally obtain
\begin{equation}
 \left({da\over dt}\right)_s = \frac{4\alpha}{9}\,{\Phi\over
 \omega_{\rm rev}}
  \,{E_{R'}\sin\delta_{R'}\over 1+\chi}\, \sin^2\gamma\; .
 \label{twentysix}
\end{equation}
This formula coincides exactly with the semimajor axis decay rate
due to the seasonal variant of the Yarkovsky effect
(see e.g. Rubincam 1998; Farinella et~al. 1998).
 
\subsection{Solution for the diurnal/mixed components}
Next, we determine the ``tesseral'' coefficients $t'_{1\pm1}(r';
\zeta)$ of the temperature development (\ref{eleven}), which are
related to the out--of--axis thermal force components $(f_X,
f_Y)$ [see (\ref{fifteen})]. Given the Fourier expansion of the
insolation coefficients $a_{1\pm 1}(\zeta)$ from (\ref{seven}) and
(\ref{nine}) we may expect
\begin{equation}
 t'_{1\pm 1}(r';\zeta) = \tau^{\pm}_{1}(r')\, \zeta^{\pm (m+1)} +
  \tau^{\pm}_{-1}(r')\, \zeta^{\pm (m-1)}\; , \label{twentyseven}
\end{equation}
and we find again from Eq. (\ref{twelve}) that the radial
amplitudes $\tau^{\pm}_{\pm 1}(r')$ satisfy a system of
decoupled homogeneous spherical Bessel equations.
Their solution reads
\begin{eqnarray}
 \tau^\pm_{1}(r') &=& c^\pm_{1}\;\; j_1\left[\sqrt{\mp i\left(m
  + 1\right)}r'\right]\; , \label{twentyeight} \\
 \tau^\pm_{-1}(r') &=& c^\pm_{-1}\, j_1\left[\sqrt{\mp i\left(m
  - 1\right)}r'\right]\; . \label{twentynine}
\end{eqnarray}
So far, we have been keeping the normalization of the radial
coordinate $r$ by the penetration depth $l_s$ of the seasonal
thermal wave. However, a normalization by the penetration depth
$l_d = l_s/\sqrt{m}$ of the diurnal wave is more suitable now,
and will be used hereafter in this Section. Thus
\begin{eqnarray}
 \tau^\pm_{1}(r') &=& c^\pm_{1}\;\; j_1\left[\sqrt{\mp i\left(1
  + 1/m\right)}r'\right]\; , \label{thirty} \\
 \tau^\pm_{-1}(r') &=& c^\pm_{-1}\, j_1\left[\sqrt{\mp i\left(1
  - 1/m\right)}r'\right] \; .\label{thirtyone}
\end{eqnarray}
After deriving the constant factors $c^\pm_{\pm 1}$ {}from the
boundary conditions we obtain
\begin{eqnarray}
 f_X+if_Y\!\! &=&\!\! -{4\over 9}\,{\alpha\Phi\over
  1+\chi} \Bigl\{\sin^2 {\gamma\over 2}\, E_{R_+'} \exp(
  -i\delta^n_{R_+'})\, \zeta^{-1} \label{thirtytwo} \\
 & &\qquad\quad\;\;\, + \cos^2 {\gamma\over 2}\, E_{R_-'} \exp(-i
  \delta^n_{R_-'})\, \zeta\Bigr\}\,\zeta^{-m}\, . \nonumber
\end{eqnarray}
Here we defined $R_\pm' = \sqrt{1\pm 1/m}\, R'$, while the remaining
quantities are the same we used previously. Note that, like
the insolation coefficients $a_{1\pm 1}(\zeta)$, the ``diurnal''
force components (\ref{thirtytwo}) also depend on a linear
combination of the rotation and revolution frequencies, but not
on the rotation frequency itself. This statement is obviously
``coordinate dependent''. To make a more obvious link
to the results of Vokrouhlick\'y (1998a; Eq.~30), one may
transform $(f_X,f_Y)$ force components to the body--centered frame
with the axis $z$ along its spin axis and the axis $x$ so that
the local direction to the Sun lies in the $xz$--plane (this
system has been used by Vokrouhlick\'y, 1998a). Denoting
the equatorial Yarkovsky force components in this reference system 
by we obtain
\begin{equation}
 f_x+if_y = {f_X+if_Y\over \sin\theta_0}
  \left(\sin^2{\gamma\over 2}\;\zeta + \cos^2{\gamma\over 2}\;
  \zeta^{-1}\right)\;\zeta^m\; . \label{thirtytwobis}
\end{equation}                  
Fourier development of (\ref{thirtytwobis}) contains
infinite series of spectral lines depending on both rotation
and revolution frequencies. In what follows, we use
(\ref{thirtytwo}) because of its simplicity.
 
Like in the previous Section, we are primarily interested in the
contribution of the $(f_X,f_Y)$ force components to the
along--track perturbation on a quasi--circular orbit. Simple
algebra yields
\begin{eqnarray}
 T_2 &=& {1\over 2} \left(f_Y-if_X\right)\, \left(\sin^2{\gamma\over
  2}\; \zeta^{m+1} - \cos^2{\gamma\over 2}\;\zeta^{m-1}\right)\, +\,
  {\rm C.C.} \nonumber \\ & & + {\cal O}(e)
 \label{thirtysix}
\end{eqnarray}
for the along--track perturbation force; C.C. means the complex
conjugate quantity of the previous term.
Taking the average of $T_2$ over one revolution around
the Sun we obtain the mean rate of change of the semimajor axis
\begin{eqnarray}
 \left({da\over dt}\right)_d &=& -{8\alpha\over 9\omega_{\rm rev}}\,{
  \Phi\over 1+\chi} \Bigl\{\cos^4 {\gamma\over 2}\, E_{R_-'}
  \sin\delta_{R_-'} \nonumber \\
 & &\qquad\qquad\qquad - \sin^4 {\gamma\over 2}\, E_{R_+'} \sin\delta_{R_+'}
  \Bigr\}\, . \label{thirtyseven}
\end{eqnarray}
where the index $d$ reminds us that we are dealing with the diurnal
components $(f_X,f_Y)$ of the thermal force. Interestingly,
Eq. (\ref{thirtyseven}) shows that the diurnal effect
on the semimajor axis is not simply proportional to the
cosine of the obliquity -- as most commonly used -- but depends
on $\gamma$ in a more complicated way. However, the classical
$\cos\gamma$ result can be recovered realizing that in
the typical astronomical applications $\omega_{\rm rot}/\omega_{\rm rev}
= m \gg 1$. Then, to a high accuracy we can set $R_\pm' \simeq R'$
and $\cos^4
\gamma/2 -\sin^4 \gamma/2 = \cos\gamma$. Then, Eq. (\ref{thirtyseven})
becomes identical with the classical result (e.g. Vokrouhlick\'y
1998a), confirming that in the fast rotation limit ($m\rightarrow
\infty$) the diurnal variant of the Yarkovsky effect is naturally
decoupled from its seasonal counterpart. It is easy to check
that in all astronomical applications listed in Section 1 the
exact result (\ref{thirtyseven}) differs from the classical formula
by less than one part in $10^3$.
 
Another consequence of Eq. (\ref{thirtyseven}) is that the
diurnal effect on the semimajor axis does not vanish
exactly at $90^\circ$ obliquity. This asymmetry implies that
the corresponding semimajor axis drift does not average out to zero
when an isotropic distribution of spin axes is assumed
(e.g., due to frequent and random impact reorientation events).
Denoting by angled brackets such an average over all the spin
orientations, we obtain
\begin{eqnarray}
 \biggl<\left({da\over dt}\right)_d\biggr>\!\! &=&\!\! {8\alpha\over 27
  \omega_{\rm rev}}{\Phi\over 1+\chi} \Bigl\{E_{R_+'}\sin\delta_{R_+'}\!\!
  - E_{R_-'} \sin\delta_{R_-'}\Bigr\}\, . \nonumber \\ & &
  \label{thirtyeight}
\end{eqnarray}
However, unless the $m$ parameter is unrealistically small,
collisional events will tend to spin up the rotation rate besides
reorienting the spin axis, and in this case the long--term
average (\ref{thirtyeight}) of the diurnal semimajor axis effect is
very small. Quantitatively, one can easily show that
$\langle (da/dt)_d \rangle \propto 1/m$ at large $m$.
 
\section{Conclusions}
The main results of this paper can be summarized as follows:
\begin{itemize}
\item We have developed a self--consistent, unified linear theory for
 the thermal Yarkovsky effects on a spherical, rotating body which
 orbits around a radiation source. The solution contains
 consistently all terms that depend on both the rotation and the
 revolution frequency and their mutual combinations. In the
 limit of a rapidly rotating body the two classical
 variants of the thermal effects -- the so--called diurnal and the
 seasonal effects previously modeled in a separate fashion -- are
 recovered at the same time.
\item Besides rederiving the previously known results for the seasonal
 and diurnal variants of the Yarkovsky effect, we have computed and
 clarified the contribution of the ``mixed'' terms to the mean
semimajor axis rate of an orbiting body.
\item We have found that in most
 relevant astronomical applications, such as the motion of asteroidal
 fragments in the main asteroid belt, the relative contribution of the
 newly derived mixed terms to the Yarkovsky semimajor axis
 perturbations is less than one part in $10^3$.  This is a negative
 result, but it is important to justify a number of recent studies
 taking into account only the limit variants of the Yarkovsky effect,
 namely the diurnal and seasonal effects, in dealing with the
 meteorite delivery issue.
 Only in the case of very slowly rotating bodies orbiting in the inner
 Solar System (a possible case would be that of Mercurian ejecta
 escaped from the planet's gravity field), the more exact solution derived
 in this paper might be used to improve the accuracy of the results.
\end{itemize}
 
\acknowledgements{The author thanks P.~Farinella for helpful discussions
 and improvements in the style of the final version of the paper
 and D.P.~Rubincam (the referee) for pointing out an error in the
 submitted manuscript.  This study has been finished at
 OCA/CERGA (France) thanks to the summer 1998 research grant and at
 the University of Pisa (Italy) thanks to a ``visiting scientist''
 grant of the ``Gruppo Nazionale di Astronomia'' of the Italian National
 Research Council (CNR).}
 
\bigskip
\noindent{\bf References}
\smallskip
\begin{description}
\item Farinella P., D. Vokrouhlick\'y, 1996, Planet. Sp. Sci. 44,
 1551
\item Farinella P., D. Vokrouhlick\'y, 1998, submitted
\item Farinella P., D. Vokrouhlick\'y, W.K. Hartmann, 1998,
 Icarus 132, 378
\item Hartmann W.K., P. Farinella, D. Vokrouhlick\'y, S.J.
 Weidenschilling, A. Morbidelli, F. Marzari, D.R. Davis,
 E. Ryan, 1998, Meteoritics Planet. Sci., in press
\item Morbidelli A., Gladman, B., 1998, Meteoritics Planet. Sci. 33,
 999
\item Rubincam D.P., 1994, Theor. Appl. Climatol. 48, 195
\item Rubincam D.P., 1995, J. Geophys. Res. 100, 1585
\item Rubincam D.P., 1998, J. Geophys. Res. 103, 1725
\item Vokrouhlick\'y D., 1998a, AA 335, 1093
\item Vokrouhlick\'y D., 1998b, AA 338, 353
\item Vokrouhlick\'y D., P. Farinella, 1998a, AA 335, 351
\item Vokrouhlick\'y D., P. Farinella, 1998b, AJ 116, 2032
\end{description}
\end{document}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
