%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Small corrections by PF, 6 February 1998
% Korekce postrehnute PO submitovani, DV Feb 27 98
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\documentstyle[epsf]{l-aa}
\newcommand{\etal}{{\it et~al.}}
\newcommand{\etall}{{\it et~al. }}
\newcommand{\boldM}[1]{\mbox{\boldmath{${#1}$}}}
%\renewcommand{\baselinestretch}{2}
%\onecolumn
\thesaurus{05.03.1; 07.13.2}
\title{Diurnal Yarkovsky effect as a source of mobility of
 meter-sized asteroidal fragments}
\subtitle{I. Linear theory}
\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 \\
    (Electronic mail: vokrouhl@mbox.cesnet.cz)}
\date{\today}
\begin{document}
\maketitle
%\offprints{D. Vokrouhlick\'y}
\markboth{D. Vokrouhlick\'y: Diurnal Yarkovsky effect ...}{}
\begin{abstract}
 A linear theory for the heat conduction in a spherical, solid and
 rotating body illuminated by solar radiation is developed in
 detail. The principal aim is to compute the recoil force, due to
 thermally reemitted radiation, which is commonly known as the
 ``Yarkovsky force''. We concentrate on the thermal effect which
 depends on the rotational period of a body rather than on the period
 of revolution around the Sun and deal with the general case of an
 arbitrary obliquity of the spin axis to the orbital plane.
 This ``diurnal'' thermal effect is considered to be an important
 source of mobility for meter-sized stony asteroidal fragments
 in the main belt. We compare our results with those of previous
 authors and show that the results of Peterson (1976) are accurate for
 meter-sized asteroidal bodies (although he used unrealistically long
 rotation periods).
\keywords
{Celestial mechanics -- Minor planets, asteroids -- Meteoroids}
\end{abstract}
 
\section{Introduction}
 
The Earth's vicinity in the solar system is by no means an empty
space.  A number of small objects of sizes ranging from kilometers
down to micrometers, with a typical power-law distribution, cross the
Earth orbit (Rabinowitz 1993, 1994; Ceplecha 1992, 1996).  It has been
recognized long ago that these objects are relatively short-lived, and
are being permanently replaced by a population of similar objects from
the outer regions of the solar system.  Their most prominent sources
are the asteroid belt and the short-period comets.  Transport
mechanisms by means of which such bodies can be delivered to the
Earth's neighbourhood have been studied intensively over the last
decades.  The crucial point has been the recognition of a special role
of the mean-motion and secular resonances in the main asteroid belt
(e.g.  Morbidelli et~al. 1994, Gladman et~al. 1997, Migliorini et~al.
1997).  However, several issues remain to be understood in order to
fit all observational data.  In particular, it has been recently
argued that the fast transport by resoant effects (on time scales of a
few Myr) must be necessarily preceded by a relatively slower phase (on
time scales of $10-100$~Myr).  The latter may be dominated by orbital
perturbations due to the thermal Yarkovsky effects (Farinella et~al.
1998, Hartmann et~al. 1998), which motive us to reconsider the
importance and develop a precise modelling of these thermal orbital
effects.
 
The present series of papers is devoted to a special type of
thermal effects -- to so called diurnal effect, which depends on the
rotational period of a body. Farinella et~al. (1998)
concluded that this effect is most relevant for the delivery of two
important classes of bodies: (i) submeter stone bodies which are
the most frequent parents of meteorites, and (ii) $20$ --
$70$~meter regolith covered bodies which, interestingly, are
overabundant among the near-Earth objects (Rabinowitz 1993,
1994). A first treatment of the problem of the orbital
perturbations due to the diurnal effect was given by Vokrouhlick\'y and
Farinella (1998a) (see also Afonso et~al. 1995, but notice discussion
below). They considered the interaction of these objects with the
$\nu_6$ resonance at the inner edge of the asteroid belt. Our
present aim is to develop a complete numerical model allowing
to include both gravitational interactions with all the planets and a
full treatment of the diurnal Yarkovsky effect in any region of
the asteroid belt. A first issue we face in this program is a
proper modelling of the instantaneous thermal force. This paper
develops such a theory. First, however, let us briefly recall the
previous work on the Yarkovsky effect in celestial mechanics.
 
The modern history of the applications of the thermal recoil force in
solar system dynamics begins with a classical paper by \"Opik
(1951).  He is also at the origin of an interesting ``mystery''
surrounding the effect.  Referring to an unpublished pamphlet of an
unknown Polish engineer Yarkovsky, \"Opik decided to name this effect
{\it Yarkovsky effect} in his honour.  Without performing any precise
calculations, \"Opik gave some basic estimates for the magnitude of
the effect.  Interestingly, if one disregards numerical factors of the
order of unity, \"Opik's results are in full agreement with today's
models, at least for the diurnal effect (see below).
 
Radzievskii (1952) then presented the first detailed, mathematical
study of the thermal effects on a sphere.  Unfortunately, his work is
limited by two unrealistic assumptions: (i) the body's spin axis must
be normal to the orbital plane, and (ii) the body is ``large'', when
compared to the depth of a thermal wave.  These simplifications made
his calculations of little use for a further progress; nevertheless
his work remains a pioneering one.
 
Radzievskii's work was followed by Peterson (1976).  Though his
analysis was much more detailed, he actually focused on aspects of the
thermal effects which have little relevance: for instance, he analyzed
in detail the complicated second-order solution of the heat conduction
equation in a cylinder.  However, from what we presently know in case
of the seasonal effect (e.g.  Vokrouhlick\'y and Farinella 1998a,b),
the non-linearity effects do not produce any new qualitative feature.
For ``large'' bodies they alter the amplitude of the thermal effects
by about $15$--$20\%$, whereas in the case of ``small' bodies the
difference between the linear and non-linear solutions is below the
$1\%$ level.  Peterson realized the importance of deriving at least
some estimates for spherical bodies, and invented a clever (though
somewhat cumbersome) method of converting the results for a cylinder
to a sphere.  His results, like those of Radzievskii, were correct for
``large'' bodies only, and he had to follow the qualitative arguments
of \"Opik (1951) to make any estimates for the ``small'' bodies.  Yet,
his work was the frist to discuss in detail the possible applications
of the Yarkovsky effect in meteoritics and solar system dynamics.  He
was almost alone in advocating the relevance of non-graviational
dynamical mechanisms at that time, and the introduction to his 1976
paper still gives an interesting historic and conceptual discussion.
Influenced by Peterson's work, Burns et~al.  (1979) reserved a special
chapter to the thermal effects in their comprehensive study of
radiative forces in the solar system.
 
Despite the progress mentioned above, one major point had been missing
in all these studies of the thermal effects (although the first
indications of new concepts may be traced back to Burns et~al. 1979).
Interestingly, a key inspiration for a generalization of the classical
Yarkovsky effect came from the investigation of the motion of the
artificial satellite LAGEOS, whose orbit is perturbed by thermal
effects similar to those affecting the small natural bodies in the
solar system (see e.g.  Rubincam 1987, Vokrouhlick\'y and Farinella
1998a).  The reason is that LAGEOS has an important well-known
property: a very fast rotation (e.g.  Farinella et~al. 1996).
Rubincam (1987) thus developed a ``LAGEOS-taylored'' technique for
computing thermal force perturbations on a rapidly spinning body.  If,
say in 1990, we had tried to compare models of thermal effects acting
on natural solar system bodies and models applied to LAGEOS, we would
have beeb surprised by a number of differenct assumptions and
techniques.  Yet, in both cases the physics is essentially the same.
 
Rubincam (1995) first tried to reconcile the models of thermal effects
in the two cases.  He applied a LAGEOS-like modelling to the case of
asteroidal fragments, after having checked the physical consistency.
A similar approach led Farinella et~al.  (1998) to reconsider the
concepts of large vs. small size and slow vs. fast rotation of the
fragments, both introduced by Burns et~al.  (1979).  The variety of
possible thermal effects was found to be larger than recognized
before.  Regarding the orbital effects, a new classification based on
``diurnal'' vs. ``seasonal'' variants was suggested.
 
Farinella et~al.  (1998) demonstrated that the diurnal Yarkovsky
effect may play a very important role in the dynamics of stony
fragments of asteroids, especially in the applications to meteoritics
(i.e., for meter-sized bodies) and in delivering Tunguska-sized cosmic
bodies (10 to 100 m across) to the Earth's vicinity.  However, their
discussion of the diurnal effect was based on Peterson's results, and
they did not try to develop a more accurate thermal force model.  The
development of such model is the primary goal of this series of
papers.
 
In the present paper we give a complete formulation of the diurnal
Yarkovsky effect, in which we suitably linearize the boundary emission
term.  Afonso et~al.  (1995) also computed explicitly the diurnal
Yarkovsky acceleration of a spherical fragment, but only in the
assumption that the fragment's spin axis is normal to the orbital
plane.  Owing to frequent collisions of fragments, resulting in
changes of the orientation of their spin axes, the very special
configuration considered by Afonso et~al. is unrealistic and thus
their formulation is incomplete.  Since the typical time scale for the
change of the spin axis orientation of meter-sized stony fragments is
a few Myr only (Farinella et~al., 1998), in order to progress in a
numerical exploration of the orbital perturbations due to thermal
effects on meter-sized objects, we need a formulation which allows for
any mutual orientation of the spin axis and the orbital plane.  The
results which we are going to present and discuss are valid for an
arbitrarily oriented spin axis.
 
\section{Theory}
\subsection{General formulation and linearization}
 
The heat conduction in a solid medium is described by the parabolic
equation (often called Fourier equation; e.g.  Landau and Lifchitz,
1986) \begin{equation} \rho C {\partial T\over \partial t} = K
\nabla^2 T \; , \label{1} \end{equation} yielding a distribution of
temperature $T$ throughout the medium at any time $t$. $\nabla^2$ is
the Laplace operator, $K$ is the thermal conductivity, $C$ the
specific heat and $\rho$ the density of the material.  In principle
all these three parameters may be temperature dependent, which would
result in a more complex form of the heat conduction equation.
However, in this paper we shall neglect these phenomena and adopt
average quantities of the physical parameters over the temperature
range involved.
 
Equation (\ref{1}) must be supplemented by an appropriate boundary
constraint at the surface of the body and by the condition that the
temperature is regular inside the body. The boundary condition is
provided by the conservation of energy and is given by
\begin{equation}
 \epsilon\sigma T^4 + K \left({\bf n}\cdot{\partial T\over \partial
  {\bf r}}\right) = \alpha\,{\cal E} \; . \label{2}
\end{equation}
The first term on the left-hand side accounts for the energy thermally
reradiated by the body (the isotropic Lambert's law is assumed),
$\epsilon$ is the emissivity and $\sigma$ the Stefan-Boltzmann
constant; the second term gives the energy conducted to the deeper
layers of the body, with ${\bf n}$ the unit vector normal to the
surface of the body.  The right-hand side of (\ref{2}) gives the
radiation energy entering the unit surface area of the body per unit
time, with $\alpha$ the absorption coefficient and ${\cal E}$ the
external radiation flux.
 
A general solution of equations (\ref{1}) and (\ref{2}) is
fairly complicated even if the geometry of the body is simple
(plane-parallel, spherical etc.). The non-linear, fourth-power
emission law in the first term of (\ref{2}) is the main source of
difficulties. A standard technique to handle this problem
is based on the assumption that the temperature throughout
the body does not differ much from some average value. Hereafter
we shall adopt this approximation and split the temperature in the
following way: $T = T_{\rm av} + \Delta T$. The mean temperature
$T_{\rm av}$ will be chosen such that $\Delta T\ll T_{\rm av}$.
If the previous condition is fulfilled, we can linearize the emission
term as $T^4 \approx T_{\rm av}^4 + 4T_{\rm av}^3 \Delta T + ...$ .
Neglecting higher-order terms in this expansion allows us
to treat the problem analytically in proper variables,
which follow from symmetries of the body. Assuming spherical
the fragments, we shall use spherical coordinates $(r,\theta,\phi)$.
The origin $r=0$ coincides with the center of the body and the
colatitude $\theta$ is measured from the fragment's spin axis.
The origin of the $\phi$ coordinate, which will be discussed later,
depends on an appropriate coordinate frame.
 
\subsection{The most suitable variables}
Before embarking on the solution of the problem a careful choice of
variables is necessary. Following the previous
work (Spencer et~al., 1989; Vokrouhlick\'y and Farinella, 1998b,c),
we adopt the following set of the non-dimensional quantities:
\begin{itemize}
 
\item the radial coordinate $r$ will be scaled by the thermal length
 $l_s$ given by
 \begin{equation}
  l_s = \sqrt{K\over \rho C \omega}\; , \label{3}
 \end{equation}
 and shall be denoted by $r' = r/l_s$, $\omega$ being the angular
 velocity of the fragment's rotation;
\item the time $t$ will be replaced by a complex variable $\zeta$
 given by
 \begin{equation}
  \zeta = \exp(i \omega t) \label{4}
 \end{equation}
 (here $i=\sqrt{-1}$ is the imaginary unit);
\item the temperature $T$ will be scaled by an auxiliary value
 $T_\star$  defined by
 \begin{equation}
  \epsilon \sigma T_\star^4 = \alpha {\cal E}_\star\; , \label{5}
 \end{equation}
 where ${\cal E}_\star$ is the solar radiation flux at the position
 of the fragment. The resulting non-dimensional variable will be denoted
 by $T' = T/T_\star$, and similarly we define $\Delta T' = \Delta
 T/T_\star$;
\item the energy source term ${\cal E}$ in the right hand side of
 (\ref{2}) will be scaled by the reference flux ${\cal E}_\star$,
 so we define ${\cal E}' = {\cal E}/{\cal E}_\star$.
\end{itemize}
 
Some comments are appropriate here.  First, the auxiliary temperature
$T_\star$ should not be confused with the mean temperature $T_{\rm
av}$.  Later on we shall see that the most natural choice for the
scaled mean temperature is $T_{\rm av}' = T_{\rm av}/ T_\star =
1/\sqrt{2}$ (see also Rubincam 1995, 1998).  By the way, we remark
that Afonso et~al.  (1994) missed this point and chose $T_{\rm
av}=T_\star$.  Secondly, the above-mentioned choice of variables
reduces the number of parameters in the transformed equations
(\ref{1}) and (\ref{2}) to a minimum set, particularly the thermal
parameter $\Theta$,
\begin{equation}
 \Theta = {\Gamma \sqrt{\omega}\over \epsilon \sigma T_\star^3}\; ,
  \label{6}
\end{equation}
where $\Gamma = \sqrt{\rho C K}$ is the thermal inertia, and the
scaled radius of the body, $R'=R/l_s$. The thermal parameter $\Theta$
defined above differs by a factor of the order of unity from
similar quantities introduced in the literature. Notably, Peterson's
$P$-variable is given by $P = (\pi^{3/4}/ \sqrt{2})\, \Theta$
(Peterson, 1976), while Farinella et~al. (1998) define
$\Theta_\omega = (\sqrt{2}/ \pi)\, \Theta$.
 
Adopting this new set of variables, the heat conduction equation
(\ref{1}) has the following form
\begin{eqnarray}
 i \zeta {\partial \over \partial \zeta}\,& &\!\!\!\!\!\!\!\!\!
  \Delta T'(r';\theta,\phi;\zeta) = \label{7} \\&=&
   {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{8}
\end{equation}
The linearized boundary condition (\ref{2}) reads
\begin{equation}
 \sqrt{2} \Delta T'+ \Theta\left({\partial \Delta T'\over \partial r'}
  \right)_{R'} = \alpha\,\Delta {\cal E}'\; ,\label{9}
\end{equation}
where the right-hand term is defined by ${\cal E}' = \frac{1}{4} +
\Delta {\cal E}'$.  Here, the first term, ${\cal E}'_{\rm av}=1/4$, is
the averaged irradiation of the fragment's surface which determines
the averaged temperature mentioned above: $\epsilon \sigma T_{\rm
av}^4 \equiv\alpha {\cal E}_{\rm av}$. The expansion of the source
term $\Delta {\cal E}'$ deserves a special care and is related to the
reference systems which will be used in the following.  Notice that
this is a major issue which was not properly solved by Afonso et~al.
(1995), and we devote the next section to its discussion.
 
\subsection{Radiation source term}
The principal reference system in which we shall solve the heat conduction
problem is rigidly rotating with the fragment. Its $z$-axis is
aligned with the unit vector ${\bf s}$ of the fragment's spin axis and
the $x$-axis is chosen so that the Sun lies in the $xz$-plane at the
instant $t=0$. The unit position vector of the Sun, ${\bf n}_0$, in
this system has the components
\begin{equation}
  {\bf n}_0 = \pmatrix{\frac{1}{2}\zeta\sin\theta_0 \cr
                       \frac{i}{2}\zeta\sin\theta_0 \cr
                       \cos\theta_0\cr} + {\rm C.C.}
 \; , \label{10}
\end{equation}
where $\theta_0$ is the solar colatitude and ${\rm C.C.}$ is a complex
conjugate quantity.  A crucial point here is the time variability of
the quantities in (\ref{10}).  In the case of the diurnal Yarkovsky
thermal effect we {\it assume} that the solar colatitude $\theta_0$ in
(\ref{10}) is constant.  This means that the thermal response of the
fragment is local or ``instantaneous'', because the spin period is
much shorter than the orbital one.  On the contrary, in the case of
the seasonal Yarkovsky thermal effect one assumes a variable solar
colatitude and averages over $\zeta$ (i.e., the rotation of the
fragment).  These are the two extreme cases of the Yarkovsky thermal
effects.  More in general one should assume that $\theta_0$ is
time-dependent and still the source term $\Delta {\cal E}'$ is not
averaged over the fragment's rotation.  In other words, one should
deal with a generic mixture of the diurnal and seasonal Yarkovsky
effects.  In this paper we adopt the deal with the diurnal case only,
and therefore we use the approximation of a constant solar colatitude
$\theta_0$.
 
We can now proceed
directly to the evaluation of the radiation source term $\Delta
{\cal E}'$ in (\ref{9}). Given a surface element with a unit normal
vector ${\bf n}(\theta,\phi) = (\sin\theta\cos\phi,\sin\theta\sin\phi,
\cos \theta)^T$, we have
\begin{eqnarray}
 {\cal E}' &=& {\bf n}(\theta,\phi).{\bf n}_0(\theta_0,\zeta)
  \qquad\;\, \mbox{\rm if $({\bf n}.{\bf n}_0)>0$,} \nonumber \\
  &=& 0 \qquad\qquad\qquad\qquad\quad \mbox{\rm otherwise.}
  \label{11}
\end{eqnarray}
Obviously, it would be too cumbersome to work analytically with such a
piecewise-defined source function. Moreover, as we shall see below, we
do not need the whole information about ${\cal E}'$ for
computing the thermal force. It is suitable to represent
(\ref{11}) in a spherical harmonics expansion
\begin{eqnarray}
 {\cal E}' &=& {\bf n}(\theta,\phi).{\bf n}_0(\theta_0,\zeta)\nonumber
  \\ & =&\sum_{n\geq0}\,\sum_{k=-n}^n a_{nk}(\theta_0,\zeta)\,
   Y_{nk}(\theta,\phi)\; , \label{12}
\end{eqnarray}
from which only the monopole ($n=0$) and dipole ($n=1$) parts will be
relevant. One can prove the relation $a_{nk}(\theta_0,\zeta)=
b_{nk}(\theta_0)\,\zeta^k$, a particular case of which is
\begin{equation}
 b_{00} = a_{00} = \frac{\sqrt{\pi}}{2} \label{13}
\end{equation}
for the monopole coefficient, and
\begin{eqnarray}
 b_{10}(\theta_0) &=& \sqrt{{\pi \over 3}}\, \cos\theta_0 \; ,
  \label{14} \\
 b_{1\pm 1}(\theta_0) &=& \mp\sqrt{{\pi \over 6}}\, \sin\theta_0
  \label{15}
\end{eqnarray}
for the dipole coefficients.  The best technique to obtain equations
(\ref{14}) and (\ref{15}) is to compute them for the special case
$\theta_0=0$ (the solar direction along the fragment's spin axis) by
simple quadratures.  Then the general case is determined by a
transformation of the right-hand side of (\ref{12}) to a new system,
with a given solar direction.  The formulae for the transformation of
the spherical functions $Y_{nk}$ using Wigner's matrixes are discussed
in the standard textbooks of Wigner (1959) and Edmonds (1974).  In the
astronomical context we refer, for instance, to \v Sidlichovsk\'y
(1983).
 
Recalling that $Y_{00} = 1/(2\sqrt{\pi})$,
we observe that the monopole term in (\ref{12}) represents exactly the
averaged irradiation of the fragment. The important
terms for determination of the thermal force, which contribute to
the $\Delta {\cal E}'$ source, are those of the dipole part.
 
\subsection{Regular solution satisfying the boundary condition}
After having discussed in detail the source term, we are ready to find
a solution of the problem (\ref{7}) with the surface condition
(\ref{9}).  First we determine a general form of the solution
which is regular throughout the body.
 
The linearity of the heat conduction equation (\ref{7}) allows us a
convenient separation of the variables.  Thanks to the spherical
system of coordinates we observe that an expansion in spherical
harmonics is the best choice.  Thus, we write
\begin{equation}
 \Delta T'(r';\theta,\phi;\zeta) = \sum_{n\geq 0}\,\sum_{k=-n}^n
  t'_{nk}(r';\zeta)\,Y_{nk}(\theta,\phi)  \; . \label{16}
\end{equation}
Considering the properties of the coefficients of the source terms
expansion (\ref{12}), we realize that $t'_{nk}(r';\zeta)=\tau'_{nk}(r')
\,\zeta^k$, where the radial functions satisfy the system of decoupled
ordinary differential equations
\begin{equation}
 \left\{{d \over d r'}\left(r'^2{d \over d r'}\right) -\left[n
  \left(n+1\right) + i k r'^2\right]\right\}\,\tau'(r') = 0 \; .
  \label{17}
\end{equation}
Their general solution, which is regular inside the body (i.e. at
$r'=0$), reads
\begin{eqnarray}
 \tau'_{nk}(r') &=& c_{n}\, r'^n\qquad\qquad\qquad\;\mbox{\rm for}
  \;\; k=0 \; ,\label{18} \\
 \tau'_{nk}(r') &=& c_{nk}\, j_n(\sqrt{-ik}r')\qquad\mbox{\rm for}\;\;
  k\neq 0\; , \label{19}
\end{eqnarray}
where $j_n(z)$ are the spherical Bessel functions of order $n$ (note
the complex argument $z$).
The constants $c_n$ and $c_{nk}$ are to be determined by fulfilling the
boundary constraint (\ref{9}).
 
Substituting (\ref{18}) and (\ref{19}) into (\ref{9}) we find
\begin{eqnarray}
 c_{n} &=& \!{b_{n0}\over \sqrt{2} R'^n} {1\over 1+n\lambda}\; ,\label{20}
   \\
 c_{nk} &=& \!{b_{nk}\over \sqrt{2}\, j_n\left(\sqrt{-ik}R'\right)}{1\over
  \left[1+\lambda\,{z\over j_n(z)}{d\over dz}j_n(z)\right]}_{z=\sqrt{-ik}
  R'}\; , \nonumber \\ & & \label{21}
\end{eqnarray}
where $\lambda \equiv\Theta/\sqrt{2}R'$. As noted above, a special
attention has to be paid to the dipole part, $n=1$. For this purpose
we introduce an auxiliary function $\psi(z)$, defined by
\begin{equation}
 1+\psi(z) = {z\over j_1(z)}\, {d\over dz}j_1(z)\; . \label{22}
\end{equation}
One can easily verify that $[\psi(z)/z]$ is nonsingular at $z=0$.
 
\subsection{Complete linear solution}
 
Putting together all the previous results we find that the general
solution for the surface temperature of an arbitrary element defined
by spherical angles $\theta$ and $\phi$ is given by
 
\begin{eqnarray}
 \Delta T'(R';\theta,\phi;\zeta)&=& {1\over \sqrt{2}\left(1+\lambda\right)}
  \biggl[b_{10}\left(\theta_0\right)\, Y_{10}\left(\theta,\phi\right)
  \label{23}
  \\ &+ & {b_{11}(\theta_0)\, \zeta\over 1+{\lambda\over 1+\lambda}\psi
  (\sqrt{-i}R')}\, Y_{11}\left(\theta,\phi\right) + {\rm C.C.}\nonumber \\
  &+& \mbox{\rm second and higher order terms}\biggr]\; . \nonumber
\end{eqnarray}
In the previous formula we have not given in explicit form the
quadrupole and higher multipole terms (although the required formulae
are given above, these terms will not play any role in computing the
resulting thermal force). It turns out to be convenient to
transform the previous solution (\ref{23}) into a new reference system.
This is a nonrotating frame, having the $Z$-axis aligned with the
fragment's spin vector as before. The $X$-axis is now
oriented so that the solar position lies always in the $XZ$-plane.
Spherical coordinates in the new, nonrotating system will be
denoted by $\vartheta$ and $\varphi$. Their relationship to the
previous coordinates
$\theta$ and $\phi$ is given by a trivial rotation: $\vartheta=\theta$
and $\varphi=\phi + \omega\, t$ (see Figure~\ref{f1}).
\begin{figure}[t]
\epsfxsize=8.5cm
%\epsfysize=6.cm
%\begin{center}
%\mbox{\epsfbox{aa13f1n.eps}}
%\end{center}
%\vspace{-1.cm}
\caption[]{Coordinate systems and variables introduced in the text.
 The $xyz$-system is rigidly rotating with the fragment, while
 the $XYZ$-system is fixed with respect to the Sun. The unit vectors
 are directed as
 follows: (i) ${\bf s}$ along the spin axis of the fragment, (ii)
 ${\bf n}_0$ to the local position of the Sun, and (iii) ${\bf
 n}$ normal to the considered surface element.
 The shaded region is not illuminated by sunlight.
 \label{f1}}
\end{figure}
 
Transforming the surface temperature distribution (\ref{23})
into the new coordinates we obtain
\begin{eqnarray}
 \Delta T'(R';\vartheta,\varphi)&=& {1\over \sqrt{2}\left(1+\lambda
   \right)}\biggl[b_{10}\left(\theta_0\right)\, Y_{10}\left(\vartheta,
   \varphi\right) \label{24}
  \\ &+ & {b_{11}(\theta_0)\over 1+{\lambda\over 1+\lambda}\psi
  (\sqrt{-i}R')}\, Y_{11}\left(\vartheta,\varphi\right) + {\rm C.C.}
   \nonumber \\ &+& \mbox{\rm second and higher order terms}\biggr]
   \; . \nonumber
\end{eqnarray}
Note a slight, but important difference between formulae (\ref{23})
and (\ref{24}): the latter does not depend on time through $\zeta$.
This is clearly due to the fact that in the source-oriented coordinate
system the surface temperature has to be stationary.
 
Next we introduce the auxiliary real functions $A(x)$, $B(x)$, $C(x)$,
$D(x)$ and $E(x)$, plus a phase $\delta(x)$, by
\begin{equation}
 {1\over 1+{\lambda\over 1+\lambda}\psi(z)} = {A(x)+iB(x)\over C(x)+
  iD(x)} = E(x)\exp\left[i\delta\left(x\right)\right]\; , \label{25}
\end{equation}
with $z=\sqrt{-i}R'$ and $x=\sqrt{2}R'$. After some algebra one
finds explicitly
\begin{eqnarray}
 A(x) &=& -\left(x+2\right)- e^x\left[\left(x-2\right)\cos x - x
  \sin x\right]\; , \label{26} \\
 B(x) &=& -x - e^x\left[x\cos x +\left(x-2\right)\sin x\right]\; ,
  \label{27} \\
 C(x) &=& A(x)+ {\lambda\over 1+\lambda}\times\label{28} \\ & &
   \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \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)+ {\lambda\over 1+\lambda}\times\label{29}\\& &
  \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \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 simple relations for $E(x)$ and $\delta(x)$.
 
\subsection{Thermal force and related quantities}
Having determined the temperature distribution on the fragment's
surface, we can easily compute the thermal recoil force.  Assuming
Lambert's isotropic thermal emission as before, one finds the
following components of the thermal force per unit of fragment mass,
projected onto the $X$, $Y$ and $Z$ axes of the nonrotating system
introduced above:
\begin{eqnarray} f_X + if_Y &=& -{4\alpha\over
9}\,\Phi\,{\sin\theta_0\over 1+\lambda}\, E_{R'}
\exp\left(-i\delta_{R'}\right)\; , \label{30}
\\ f_Z &=& -{4\alpha\over 9}\,\Phi\,{\cos\theta_0\over 1+\lambda}\; ,
\label{31} \end{eqnarray}
where $E_{R'} \equiv E(\sqrt{2}R')$ and $\delta_{R'} = \delta(\sqrt{2}
R')$.  The factor $\Phi \equiv (\pi R^2 {\cal E}_\star/mc)$, where $m$
is the mass of the fragment and $c$ the velocity of light, is the
usual radiation pressure parameter on a spherical body (see e.g.
Burns et~al. 1979).  We note that the remaining factors in (\ref{30})
and (\ref{31}) yield the efficiency of the thermal recoil force when
the latter is compared to the direct solar radiation pressure.
 
The physical meaning of the thermal parameter $\Theta$ has been
discussed
by Farinella et~al. (1998), who interpreted it as the ratio of
the thermal relaxation time, required for the accumulation of the
absorbed energy and its reemission, to the rotation period. As a
consequence, the limit $\Theta \rightarrow 0$, and therefore
$\lambda
\rightarrow 0$, corresponds to the case of instantaneous reemission
of the thermal energy. This is normally interpreted as a simple
diffusion of light on a macroscopic sphere (no diffraction effects;
see, for instance, Milani et~al., 1987; Vokrouhlick\'y et~al. 1993).
We can easily check that $E_{R'}\cos\delta_{R'}\rightarrow 1$ when
$\lambda\rightarrow 0$, which results into
\begin{equation}
 {\bf f}(\lambda=0) = -{4\alpha\over 9}\,\Phi\,{\bf n}_0\; ,
  \label{31bis}
\end{equation}
and obviously $f_Y(\lambda=0)=0$. The force (\ref{31bis}) is just
opposite to the local solar direction, as might have been expected,
and its magnitude is just $4/9$ of the direct radiation pressure.
\begin{figure}[t]
\epsfxsize=8.7cm
%\epsfxsize=7.5cm
%\begin{center}
%\mbox{\epsfbox{aa13f34.eps}}
%\end{center}
\vspace{-2.2cm}
 \caption[]{Coefficients $\kappa_1$, $\kappa_2$ and $\kappa_3$ vs. the
 scaled radius $R'$ of the body.  The values derived by Peterson
 (1976), given by the dashed lines, match reasonably well our results
 in the limit of ``large'' bodies ($R' \rightarrow \infty$).
 \label{f34}}
\end{figure}
 
One can see that the equatorial thermal force components
$(f_X,f_Y)$ in (\ref{30}) are rational functions of $\lambda$
(and consequently of the thermal parameter $\Theta$). A little
algebra yields
\begin{eqnarray}
 f_X = -{4\alpha\over 9}\,\Phi\,\sin\theta_0\,{1+\kappa_1\Theta\over
  1+2\kappa_1\Theta+\kappa_2\Theta^2}\; , \label{pet1}\\
 f_Y = -{4\alpha\over 9}\,\Phi\,\sin\theta_0\,{\kappa_3\Theta\over
  1+2\kappa_1\Theta+\kappa_2\Theta^2}\; , \label{pet2}
\end{eqnarray}
where the coefficients $\kappa_1$, $\kappa_2$ and $\kappa_3$ are
functions of $R'$.  Note that this result is in agreement with
Peterson's (1976) Pad\'e approximant representation [Eq.~(26) in his
paper].  However, Peterson obtained the values for these coefficients
numerically, in the limit of large bodies ($R' \approx \infty$) only.
One can easily check that there is a rather good agreement between our
results and those of Peterson, by taking the corresponding limit in
our formula (\ref{pet2}).  Figure~\ref{f34} shows the dependence on
$R'$ of the $\kappa$ coefficient (solid curves) together with the
limit values at large $R'$ obtained by Peterson (dashed lines).  Note
also a significant increase of the quadratic coefficient $\kappa_2$ in
the denominator at small $R'$, as well as a significant decrease of
the numerator coefficient $\kappa_3$ for the transverse force
component.  Both these results have to do with the important fact that
thermal effects are inhibited for small bodies, due to more efficient
heat conduction across them.  In other words, the temperature
differences in the body decay very rapidly with decreasing sizes.
 
Figures~\ref{f2} and \ref{f3} show the dependence of the
``equatorial'' force components $f_X$ and $f_Y$ on the thermal
parameter $\Theta$.  The former figure has been derived for a large
body ($R'=500$), which according to Farinella et~al.  (1998)
corresponds to a stony fragment about $1.5$~m in radius.  We note
again a good agreement with Peterson's (1976) approximation shown by
dashed lines.  As the body becomes smaller (Figure~\ref{f3},
corresponding to $R'=1$), we can remark a significant difference
between the results of the two approaches for the $f_Y$ component.
This is consistent with the previous remark concerning the behaviour
of $\kappa_2$ and $\kappa_3$ for small $R'$.
\begin{figure}[t]
\epsfxsize=8.7cm
%\epsfxsize=7.5cm
%\begin{center}
%\mbox{\epsfbox{aa13f2.eps}}
%\end{center}
\vspace{-2.2cm}
\caption[]{Normalized equatorial components of the thermal force vs.
the thermal parameter $\Theta$.  Curves 1 and 2 show
$9f_X/(4\alpha\Phi)$ and $9f_Y/(4\alpha\Phi)$, respectively.  The
radius of the body is $R'=500$, corresponding to about $1.5$~m for a
stony fragment.  Peterson's (1976) solution (dashed curves) is shown
for comparison.
 \label{f2}}
\end{figure}
\begin{figure}[t]
\epsfxsize=8.7cm
%\epsfxsize=7.5cm
%\begin{center}
%\mbox{\epsfbox{aa13f3.eps}}
%\end{center}
\vspace{-2.2cm}
\caption[]{The same as Fig.~3, but for a $R'=1$ body, corresponding to
about $10$~$\mu$m for a basalt particle.  A significant discrepancy
with respect to Peterson's results is apparent now, mainly in the
transverse force component $f_Y$, as discussed in the text.
 \label{f3}}
\end{figure}
 
Having computed the force components we can also estimate some of the
orbital effects resulting from the thermal perturbations.  In
particular the long-term variations of the semimajor axis are
essential because of the original motivation of this work: the
relevance of the Yarkovsky effect as a source of mobility for
fragments in the asteroidal belt.  To the zero order in the
eccentricity $e$, we easily obtain
\begin{equation}
 {da\over dt} = -\frac{8\alpha}{9}\,{\Phi\over n}\,{E_{R'}\sin
  \delta_{R'}\over 1+\lambda}\, \cos\gamma +{\cal O}\left(e
  \right)\; . \label{32}
\end{equation}
Note the obliquity dependent factor $\cos\gamma = {\bf s}\cdot {\bf
N}$, where ${\bf N}$ is the unit vector normal to the mean orbit; $n$
is the fragment's orbital mean motion around the Sun. We can also
remark that the $Z$-component of the thermal force (\ref{31}) does not
contribute to the formula (\ref{32}) but, obviously, yields
short-periodic perturbations of the fragment's motion.  The spin axis
of the fragment il likely to undergo random changes due to collisions
with other fragments, and if these are frequent enough we conclude
that the long-term change of the semimajor axis vanishes because
$\langle \cos\gamma\rangle =0$ (as already conluded by Burns et~al.
1979 and Farinella et~al. 1998).
 
Similarly, we can obtain the long term perturbation of the
inclination
\begin{eqnarray}
 {dI\over dt} &=& \frac{2\alpha}{9}\,{\Phi\over na}\,{s_P\left(
  1-E_{R'}\cos\delta_{R'}\right)\cos\gamma + s_Q E_{R'}\sin
  \delta_{R'}\over 1+\lambda}\nonumber \\ & & +\,{\cal O}\left(
  e\right)\; ,\label{33}
\end{eqnarray}
where $s_P = {\bf s}\cdot{\bf P}$ and $s_Q = {\bf s}\cdot{\bf Q}$ are
the projections of the spin vector ${\bf s}$ onto the position vector
${\bf P}$ of the mean pericenter and the vector ${\bf Q} = {\bf
N}\times{\bf P}$. Unlike the case of the semimajor axis, the $Z$
thermal force component (\ref{31}) does contribute to the long term
perturbation of the inclination.  Taking again the mean value of
(\ref{33}) over random orientations of the spin axis ${\bf s}$, we
conclude that the diurnal Yarkovsky effect leaves the mean inclination
unchanged over a long span of time.
 
Finally, we observe that the eccentricity's long term perturbation is
proportional to the first power of the eccentricity itself: $(de/dt)
\propto e$.  Although we might compute this term explicitly, we leave
a discussion of the analytically determined mean orbital effects to a
subsequent study.  Indeed, the main objective of this paper was the
derivation of the instantaneous thermal force components (\ref{30})
and (\ref{31}).
 
\subsection{Comparison with other approaches and discussion}
First, let us comment on the work of Afonso et~al.  (1995) whose
results have been applied to assess the diurnal thermal effects on the
orbits of spherical bodies.  It may appear that our solution is only a
slight generalization for the case of an arbitrary orientation of the
fragment's spin axis.  However, this generalization is essential for
carrying out further numerical work.  Moreover, we can see that the
Afonso et~al. solution is not really correct.  For instance, these
authors did not notice that the axially symmetric ($m=0$) dipole mode
of the temperature distribution is not time-dependent.  As a result,
the last term in their formula (31) is not correct, as it should not
depend on time.  Secondly, they neglected the difference between the
averaged temperature ($T_{\rm av}$) and the subsolar temperature
($T_\star$).  Thus a factor of two is missing in their solution.  This
fact can be easily verified in the limit of the instantaneous
rediffusion of sunlight by the fragment's surface.  Their formula (37)
yields a factor $2/9$, instead of a correct value $4/9$ derived here
[see Eq.~(\ref{31bis}) above].  As a consequence, their discussion of
the reasons for a factor two difference with respect to Peterson's
solution (see p.~790) is also incorrect.
 
\begin{figure}[t]
\epsfxsize=8.7cm
%\epsfxsize=7.5cm
%\begin{center}
%\mbox{\epsfbox{aa13f4.eps}}
%\end{center}
\vspace{-2.2cm}
\caption[]{Ratio of the maximum semimajor axis drift rate $da/dt$
(for $\gamma=0$, circular orbits) obtained by our Eq.~(\ref{32}) to
that of Peterson (\ref{34}), as a function of the radius of the body.
Thermal properties appropriate for stony fragments are assumed.  A
fairly good agreement throughout the relevant size range is observed
despite the entirely different approaches.
\label{f4}}
\end{figure}
 
A comparison with the work of Peterson (1976) is more interesting
(note that Farinella et~al. 1998 used
these results for estimating the orbital drift rate due to the
diurnal effect). Peterson's formula for the mean semimajor axis
rate, when rewritten in our variables, reads
\begin{equation}
 {da\over dt} = {\alpha\Phi\over n}\, f_P(\Theta)\, \cos\gamma
  \; ,\label{34}
\end{equation}
with the function $f_P$ given by
\begin{equation}
 f_P(\Theta) = 0.4\,{\Theta\over 1+ 0.914\,\Theta + 0.413\,
  \Theta^2}\; . \label{35}
\end{equation}
Interestingly, Peterson's result (\ref{34}) matches quite well our
general formula (\ref{32}), at least in the important case of
meter-sized stony fragments --- see Figure~\ref{f4} --- despite the
fact that Peterson used an entirely different mathematical technique.
This agrees with our previous finding -- that Peterson's results are
fairly accurate in the case of large bodies.  Note that the thermal
penetration depth corresponding to the rotation frequency is about
$l_s \approx 2.5\, \sqrt{R}$~mm if $R$ is given in meters (this
formula follows from the assumption the $\omega \propto R^{-1}$,
according to Farinella et~al. 1998).  The depth $l_s$ is even smaller
if the surface is covered by a thin layer of regolith or any porous
material which significantly decreases the thermal conductivity $K$
(Rubincam 1995; Farinella et~al. 1998).  Therefore asteroid fragments
larger than a centimeter can be safely considered ``large'', from the
point of view of the diurnal effect.
 
\section{Conclusions}
The main results of this paper may be summarized as follows:
\begin{itemize}
 
\item We have discussed in detail a linear theory for the diurnal
 Yarkovsky thermal effect on spherical rotating bodies (presumably
 asteroidal fragments). Although this is not the first paper on this
 issue, we believe it gives, for the first time, a clear and
 unambiguous solution free from the limitations and errors
 of previous studies.
\item The main parameters of the theory are: (i) the radius of the
 fragment $R$ scaled by the penetration depth of the
 thermal wave $l_s$ ($R' = R/l_s$), and (ii) the ratio between the
 thermal parameter $\Theta$ and the scaled radius $R'$ of the body.
 In the limit of large bodies, $R' \rightarrow \infty$, and the
 present solution  matches the results of Peterson (1976).
\item Given these two parameters, we derive instantaneous thermal
 force acting on the body for an arbitrary orientation of the
 spin axis.
\end{itemize}
 
Although the analytic estimates of Sec.~2.6 are given for
circular orbits, we remark that our solution (\ref{30}) and
(\ref{31}) for the thermal force is by no means limited to
low-eccentricity orbits. This is due to the simpler character of the
diurnal Yarkovsky effect, when compared with the seasonal effect.
As explained above, the diurnal effect is purely local, since at a
given
time it depends on the instantaneous dynamical variables only. On the
other hand, the seasonal Yarkovsky force depends non-locally on the
dynamical state at a number of previous instants. This complicates
the analytic expansions for eccentric orbits (see, e.g., Vokrouhlick\'y
and Farinella 1998c). In our case, the long-term analytic estimates
can be generalized to the case of arbitrary eccentricities
provided one pays attention to: (i) the expansion of the true
anomaly terms in $\sin\theta_0$ and $\cos\theta_0$, (ii) the simple
luminosity decrease relationship ${\cal E}_\star \propto (a/r)^2$ in
$\Phi$, and (iii) the orbital dependence of the thermal parameter
$\Theta$ given by
\begin{equation}
 \Theta = \Theta_0 a^{3/2} \left({a\over r}\right)^{3/2}\, .
 \label{36}
\end{equation}
Here, $\Theta_0$ is the thermal parameter at $1$~AU from the Sun,
whereas $a$ is the semimajor axis in AU and $r$ the distance from the
Sun. We leave this task for a subsequent study.
 
As far as theoretical modelling is concerned, we note that a careful
understanding of the mixture between the diurnal and seasonal
Yarkovsky effects remains a major challenge for the future.  In the
optics of this paper, one should assume that the solar colatitude
$\theta_0$ in (\ref{12}) -- (\ref{15}) is not constant but changes on
the timescale of one revolution around the Sun. The final Fourier
expansion in the time-like variable $\zeta$ would be, however, much
more complicated than in the case discussed here.
 
In the next paper, we intend to discuss some implications
of the diurnal Yarkovsky effect for meteoritics and celestial mechanics,
by applying the present results to a direct numerical integration
of the orbits of asteroid fragments.
\medskip
 
\noindent{\bf Acknowledgements}\hskip 2mm The author thanks M.~Bro\v z
for checking some calculations, J.~Bi\v c\'ak and P.~Farinella for
discussions and language corrections.  Partial support from the Czech
Grant Agency under contract No.~205/96/K119 is also acknowledged.
 
\bigskip
\noindent{\bf References}
\smallskip
\begin{description}
\item Afonso G., R.S. Gomes, M.A. Florczak, 1995, Planet. Space
 Sci. 43, 787
\item Burns J.A., P.L. Lamy, S. Soter, 1979, Icarus 40, 1
\item Ceplecha Z., 1992, AA 263, 361
\item Ceplecha Z., 1996, AA 311, 329
\item Edmonds A.R., 1974, Angular momentum in quantum mechanics, Princeton
 Univ. Press, Princeton
\item Farinella P., D. Vokrouhlick\'y, F. Barlier, 1996,
 J. Geophys. Res. 101, 17861
\item Farinella P., D. Vokrouhlick\'y, W.K. Hartmann, 1998,
 Icarus, in press
\item Gladman B.J., F. Migliorini, A. Morbidelli, V. Zappal\`a,  P.
 Michel, A. Cellino, Ch. Froeschl\'e, H.F. Levison, M. Bailey, M.
 Duncan, 1997, Science 277, 197
\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., submitted
\item Landau L.D., E. Lifchitz, 1986, {\it Hydrodynamics}, Nauka,
 Moscow (in russian)
\item Migliorini F., A. Morbidelli, V. Zappal\'a, B.J. Gladman, M.E.
 Bailey, A. Cellino, 1997, Meteoritics Planet. Sci. 32, 903
\item Milani A., A.M. Nobili, P. Farinella, 1987, Non-gravitational
 Perturbations and Satellite Geodesy, A. Hilger, Bristol
\item Morbidelli A., R. Gonczi, Ch. Froeschl\'e, P. Farinella, 1994,
 AA 282, 955
\item Peterson C., 1976, Icarus 29, 91
\item \"Opik E. J., 1951, Proc. Roy. Irish Acad. 54, 165
\item Rabinowitz D.L., 1993, ApJ 407, 412
\item Rabinowitz D.L., 1994, Icarus 111, 364
\item Radzievskii V.V., 1952, Astron. Zh. 29, 162 (in Russian)
\item Rubincam D.P., 1987, J. Geophys. Res. 92, 1287
\item Rubincam D.P., 1995, J. Geophys. Res. 100, 1585
\item Rubincam D.P., 1998, J. Geophys. Res. 103, 1725
\item Spencer J.R., L.A. Lebofsky, M.V. Sykes, 1989, Icarus 78, 337
\item \v Sidlichovsk\'y M., 1983, Bull. Astron. Inst. Czechosl. 34,
 65
\item Vokrouhlick\'y D., P. Farinella, 1997, Planet. Space Sci.
  45, 419
\item Vokrouhlick\'y D., P. Farinella, 1998a, AA, in press.
\item Vokrouhlick\'y D., P. Farinella, 1998b, Astron. J., submitted.
\item Vokrouhlick\'y D., P. Farinella, 1998c, Icarus, submitted.
\item Vokrouhlick\'y D., P. Farinella, D. Lucchesi, 1993,
 Celest. Mech. 57, 225
\item Wigner E.P., 1959, Group Theory, Academic Press, New York
\end{description}
\end{document}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

