%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% (date: Mar 2, 1998)
% English improved by PF, April 9th 1998
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\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{II. Non--sphericity effects}
\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 ..., II.}{}
\begin{abstract}
The diurnal Yarkovsky effect, a perturbing force on meter--sized
orbiting bodies caused by anisotropic heating and emission of thermal
radiation from small spinning solar system bodies, is computed for
spheroidal--shaped bodies with an arbitrary flattening.  The solution
is derived in analytical form for a body which rotates around its axis
of symmetry.  A numerical solution is presented for the more general
case of a precessing body, where the symmetry axis tumbles around a
fixed angular momentum vector.  In both cases, the obtained Yarkovsky
force is compared with the corresponding thermal force acting on a
spherical body of the same mass.  Differences of up to a factor of $3$
are found, depending on the geometric parameters of the body and its
rotational state.  However, the agreement between the derived force
values in the two cases is much better when an average over all
possible orientations of the spin axis and the solar position is
considered.  This result suggests that the simplified formulation of
the diurnal Yarkovsky force based on a spherical geometry can be used,
without introducing major mistakes, for long--term integrations of the
orbits of small asteroidal fragments in the Solar System.
\keywords
{Celestial mechanics -- Minor planets, asteroids -- Meteoroids}
\end{abstract}
 
\section{Introduction}
Thermal forces due to anisotropic heating and re-emission of thermal
radiation from small, spinning Solar System bodies (usually called
``Yarkovsky effects'') have been recently found to provide an
efficient transport mechanism for meter--sized asteroidal fragments
into resonant channels which deliver them to the near--Earth space
(Farinella et al. 1998, Hartmann et al. 1998).  Although so far most
models for these effects have assumed spherically shaped bodies, real
asteroidal fragments probably have very irregular shapes.  There are
several arguments for this conclusion: (i) fragments from asteroidal
catastrophic collisions probably have shapes similar to fragments
obtained in laboratory hypervelocity impact experiments, which
typically show approximatelly $2:1$ ratio between the longest and shortest
axes (Fujiwara et al. 1978, 1989; Catullo et al. 1984; Giblin et al.
1994); (ii) asteroid lightcurve photometry, radar observations
and spacecraft reconnaissance, all show that asteroids in the size
range between $\approx 1$~km and a few hundred km have irregular
fragment--like shapes (see e.g.  Binzel et al. 1989, Ostro et~al.
1995; Zappal\`a \& Farinella 1997); (iii) meteorites recovered on the
ground have irregular shapes, including the largest ones (a few m
across), such as Hoba and Willamette.  Note that the size range for
which these data are relevant (from $\approx 1$ cm to $>1$~km)
brackets the size range for which Yarkovsky effects are significant
(from about $0.1$ to $100$~m, see Farinella et al. 1998).
 
All of this casts doubt on the current thermal models from which
Yarkovsky effects are estimated, as they are based on questionable
simplifying assumptions as far as the shape of the body is concerned.
Actually, complete models including thermal conduction across the body
have been developed only for spheres (Vokrouhlick\'y 1998 -- hereafter
also called paper~I; see also
Afonso et~al. 1995 and Rubincam 1998).  The reason is that the view
has prevailed that a generalization for different shapes would require
a much more involved mathematical technique and the results would be
algebraically very complicated.
 
Luckily, there is an important circumstance which helps reducing the
complexity of models for the diurnal Yarkovsky effect.  In most
relevant situations the thermal penetration depth, characterizing
approximately the thickness of the surface layer undergoing
temperature variations due to external heating, is very small.  For
instance, thermal parameters appropriate for a non--porous basalt
surface would correspond to a thermal penetration depth of
about $0.4\sqrt{P}$~mm (here $P$ is the spin period of the body in
seconds), whereas an insulating regolith cover would further decrease
this figure (Farinella et~al. 1998, Hartmann et~al. 1998).  This
implies that the thermal conduction across the body is extremely weak
and each surface element can be seen as having its own thermal history
disconnected from the other elements.  In these conditions, thermal
conduction can be modelled as a one--dimensional problem for each
individual surface element.  Clearly, the body's shape does not play
an important role in this formulation.  Taking into account only the
thermal history of individual surface elements, as outlined above, one
may thus succeed in modelling Yarkovsky thermal effects for an
arbitrary shape.  Of course, when the shape has some symmetry the
required algebra becomes much less complicated.
 
It should be noted that this approach based on individual surface
elements was already used by Radzievskii in his pioneering 1952 paper
about the thermal force on a spherical body (Radzievskii 1952).  Later
on, Peterson (1976) adopted a similar formulation when developing a
theory for the thermal force on a cylinder and the same was done by
Rubincam (1995) in modeling the seasonal Yarkovsky force on a sphere.
Vokrouhlick\'y \& Farinella (1998) included the nonlinear
re--radiation effects by developing a suitable numerical model.
However, apart from the papers of Peterson (1976) and Vokrouhlick\'y
\& Farinella (1998), in all the other studies the body was assumed to
be spherical.  Thus, relatively little is known about the role of
non--sphericity effects in the Yarkovsky force theory.
 
The aim of this paper is to explore the role of the body's shape on
the resulting diurnal Yarkovsky force.  In particular, we shall
investigate how much the exact results for non--spherical bodies are
mismodelled when one adopts a simplified formulation assuming a
spherical geometry.  Such formulations are common in the studies of
the long--term dynamics of the asteroidal fragments (e.g. Afonso
et~al. 1995; Bottke et~al. 1998).  In order to assess the validity
of the corresponding
results, we need to understand the limitations of the thermal models
from which the Yarkovsky force is derived.
 
In order to get a first guess on the problem, we may use Peterson's
(1976) results.  Peterson computed the thermal force components
$(f_X,f_Y)$ perpendicular to the axis of a cylinder.  Rewriting them
into the notations introduced in this series of papers we obtain
\begin{equation}
 f_X + if_Y = -{\pi\over 3}\,\alpha\,{{\cal E}_\star Rl\over mc}\,{\sin
  \theta_0}\, {(1+\lambda)+i\lambda\over 1+2\lambda+2\lambda^2}
 \; , \label{1}
\end{equation}
with
\begin{equation}
 \lambda = {\Theta\over 4\sqrt{2}}\,\left({\pi\over \sin\theta_0}
  \right)^{3/4}  \; . \label{2}
\end{equation}
The parameters and constants are as follows: $R$ is the radius of the
cylinder, $l$ its length, $m$ its mass, ${\cal E}_\star$ the external
radiation flux, $c$ the velocity of light and $\alpha$ the optical
absorption coefficient of the surface.  The thermal parameter
$\Theta$, to be defined below, is used here instead of Peterson's
parameter $P$.  We have slightly generalized Peterson's result
[Eqs.~(17) in his paper], allowing for an arbitrary angle $\theta_0$
between the cylinder's axis and the direction toward the Sun (Peterson
assumed $\theta_0=90^\circ$).  Assuming a value $\gg 1$ for the
thermal parameter $\Theta$, we can compare (\ref{1}) to the
corresponding thermal force components for a spherical body
(\ref{23}).  Considering a cylinder of radius equal to the radius of
the sphere and a length such that the two bodies have equal mass, we
obtain that both components from (\ref{1}) are a factor
$2\sqrt{2}(\sin\theta_0/\pi)^{3/4}$ larger than the corresponding
components from Eq.~(\ref{23}) below.  When $\theta_0 = 90^\circ$,
with the solar direction perpendicular to the spin axis of the body,
this means that the force acting on a cylinfer is about $20$~\%
greater than that on a sphere with the same mass.  When $\theta_0$
becomes smaller, the solar direction gets aligned with the cylinder's
axis, and both equatorial thermal force components become much smaller
than for the comparison sphere.
 
Is this result representative enough?  Is the ``equatorial'' $\approx
20$\% error estimate valid for bodies with different shapes?  These
are the questions which will be addressed in this paper.  Of course,
we will study a specific class of on--spherical bodies, namely a
sequence of axisymmetric ellipsoids (spheroids) with an arbitrary
flattening.  Both oblate and prolate spheroids will be considered.
For the sake of simplicity, we shall first assume (Sec.~2.2) that the
rotation axis is aligned with the axis of symmetry.  In this
particular case, the solution can be obtained in a closed analytical
form.  However, real asteroidal fragments may survive for relatively
long times in a tumbling rotational state, characterized by a
misalignment between the rotation axis and the symmetry axis (e.g.,
Burns \& Safronov 1973, Harris 1994, Giblin \& Farinella 1997).  For
this reason, we will also compute the diurnal Yarkovsky force on a
body with an arbitrarily oriented spin axis.  Due to the algebraic
complexity, these latter results will be given numerically.  A linear
approximation for the grey--body thermal re-emission will be assumed
throughout this paper (see Vokrouhlick\'y 1998).
 
\section{Theory}
\subsection{Spherical bodies}
We start our discussion of the diurnal Yarkovsky effect by recalling
the results for a spherical body.  Although this case was dealt with
in paper~I, a summary is useful for two different reasons: (i) we
intend to introduce the computation method used later on in this paper
by illustrating it in the algebraically simpler case of a spherical
body, and (ii) the corresponding results will be used for a comparison
with the more general results to be derived later.  Interestingly, we
shall also demonstrate that our method overcomes an intricate and
laborious feature of Peterson's (1976) approach, namely the use the
Pad\'e approximation for his infinite series (25).  By a more direct
and simpler computation we will obtain a more precise result, because
the Pad\'e approximant method does not yield the required coefficients
with a very high precision.
 
\subsubsection{Formulation of the problem}
The one--dimensional heat conduction process in a solid medium is
described by the parabolic equation
\begin{equation}
 \rho C {\partial T\over \partial t} = K {\partial^2 T\over \partial x^2}
 \; , \label{3}
\end{equation}
yielding the distribution of temperature $T$ at depth $x$ and time
$t$.  We assume that $x$ increases toward the centre and is $0$ on the
surface. $K$ is the thermal conductivity, $C$ the specific heat and
$\rho$ the density of the material.  As in paper~I, we shall assume
that all these parameters are constant.
 
The solution of (\ref{3}) is constrained by a regularity condition
at ``infinite depth'' in the body ($x\rightarrow \infty$), and
by the energy conservation equation
\begin{equation}
 \epsilon\sigma T^4 - K \left[{\partial T\over \partial x}\right]_0
  = \alpha\,{\cal E} \; , \label{4}
\end{equation}
at the surface ($x=0$).  We use a grey--body approximation, with an
emissivity $\epsilon$ for the body's thermal re--radiation of the
absorbed energy; $\alpha$ is the absorption coefficient, ${\cal E}$
the external radiation flux on the chosen surface element and $\sigma$
the Stefan--Boltzmann constant.  The second term on the left--hand
side corresponds to the thermal energy conducted into the body.
 
The linearization of the emission term in (\ref{4}) is necessary to
obtain a sufficiently simple analytical result.  In addition to the
assumption about the ``large'' size of the body discussed above, this
linearization provides the second fundamental simplification of our
method.  Mathematically, it means splitting the temperature $T$ into a
properly chosen mean value $T_a$ and a variation $\Delta T$: $T = T_a
+ \Delta T$.
 
\subsubsection{A suitable choice of variables}
Before working out the solution of the problem, a careful choice of
variables is necessary.  Following previous work in this field
(Spencer et~al. 1989, Vokrouhlick\'y \& Farinella 1998, Vokrouhlick\'y
1998), we adopt the following set of non--dimensional quantities:
\begin{itemize}
\item the linear coordinate $x$ will be scaled by the diurnal thermal
 length $l_s$ given by
 \begin{equation}
  l_s = \sqrt{K\over \rho C \omega}\; , \label{5}
 \end{equation}
where $\omega$ is the angular velocity of rotation,
and will be denoted by $x' = x/l_s$,
\item the time $t$ will be replaced by the complex variable $\zeta$
given by
 \begin{equation}
  \zeta = \exp(i \omega t) \label{6}
 \end{equation}
 ($i=\sqrt{-1}$ is the imaginary unit);
\item the temperature $T$ will be scaled by an auxiliary parameter
$T_\star$ defined by
 \begin{equation}
  \epsilon \sigma T_\star^4 = \alpha {\cal E}_\star\; , \label{7}
 \end{equation}
 where ${\cal E}_\star$ is the solar radiation flux at the distance
 of the fragment.
 The resulting non--dimensional variable will be denoted
 by $T' = T/T_\star$, and similarly $\Delta T' = \Delta T/T_\star$;
\item the energy source term ${\cal E}$ in the right--hand side of
(\ref{4})
 will be scaled by the reference flux ${\cal E}_\star$, thus defining
 ${\cal E}' = {\cal E}/{\cal E}_\star$.
\end{itemize}
 
Before formulating the basic equations in the new set of variables,
some comments are in order.  First, the mean temperature $T_a$ defined
previosuly is not to be confused with the auxiliary value $T_\star$.
On the contrary, we shall note below that for a spherical body a
natural choice is $T'_a = T_a/T_\star = 1/ \sqrt{2}$.  Second, a major
difference between our method and that of Peterson (1976) is our
choice of a single mean temperature for {\it all} the surface elements
on the sphere.  Peterson, following the original method of Radzievskii
(1952), decided to define the mean temperature separately at each
latitude on the sphere.  However, we shall demonstrate below that this
choice turns out to be algebraically very laborious and unnecessary.
Third, the only parameter appearing in our transformed equations will
be:
\begin{equation}
 \Theta = {\Gamma \sqrt{\omega}\over \epsilon \sigma T_\star^3}\; ,
  \label{8}
\end{equation}
where $\Gamma = \sqrt{\rho C K}$ is the thermal inertia.  In agreement
with usual terminology we shall call $\Theta$ just ``thermal
parameter''.  Note that Peterson's parameter $P$ is related to
$\Theta$ by a numerical factor of order unity (e.g.  Vokrouhlick\'y
1998).
 
Adopting this new set of variables, the heat conduction equation (\ref{3})
takes the following form
\begin{equation}
 i \zeta {\partial \over \partial \zeta}\, \Delta T'(x';\theta,\phi;\zeta)
  = {\partial^2 \over \partial x'^2}\,\Delta T'(x';\theta,\phi;\zeta)
  \label{9}
\end{equation}
whereas the linearized boundary constraint (\ref{4}) reads
\begin{equation}
 \sqrt{2} \Delta T'- \Theta\left[{\partial \Delta T'\over \partial x'}
  \right]_0 = \alpha\,\Delta {\cal E}'\; .\label{10}
\end{equation}
The energy--source function $\Delta {\cal E}'$ is defined by: $\Delta
{\cal E}'={\cal E}'-\frac{1}{4}$.  The $1/4$ term here compensates for
the definition of the mean temperature, $T_a'=1/ \sqrt{2}$.
 
\subsubsection{Radiation source term}
Before we derive a solution of the system (\ref{9}) and (\ref{10}),
we comment on a suitable expansion of the source term $\Delta{\cal E}'$.
This issue is related to a definition of the reference system to
be used hereafter.
 
Considering a spherical shape of the body, we shall parametrize the
individual surface elements by the spherical coordinates $\theta$ and
$\phi$ (the colatitude $\theta$ is measured from the body's spin
axis).  The reference system used for the solution of the heat
conduction problem is rigidly rotating with the body.  Its $z$-axis is
aligned with the spin axis and the equatorial $x$-axis is chosen so
that the direction toward the radiation source (the Sun) lies in the
$xz$-plane at an arbitrary time $t=0$ (for a diagram see Fig.~1 in
Vokrouhlick\'y 1998).  Then, the unit position vector of the Sun reads
\begin{equation}
  {\bf n}_0 = \pmatrix{\frac{1}{2}\sin\theta_0 \zeta \cr
                       \frac{i}{2}\sin\theta_0 \zeta \cr
                       \cos\theta_0\cr} + {\rm C.C.}
 \; , \label{11}
\end{equation}
where $\theta_0$ is the solar colatitude and C.C. stands for a complex
conjugate quantity. The scaled solar radiation flux onto a
chosen surface element with a position unit vector ${\bf n}(\theta,
\phi)$ is
\begin{eqnarray}
 {\cal E}' &=& {\bf n}(\theta,\phi)\cdot{\bf n}_0(\theta_0,\zeta)
  \qquad\;\, \mbox{\rm if $({\bf n}\cdot{\bf n}_0)>0$,} \nonumber \\
  &=& 0 \qquad\qquad\qquad\qquad\quad \mbox{\rm otherwise.}
  \label{12}
\end{eqnarray}
This function can be expanded into a Fourier series
\begin{equation}
 {\cal E}' = {\sum_{n\in{\bf Z}}}\, \varepsilon_n (\theta,\phi;
  \theta_0)\,\zeta^n\; , \label{13}
\end{equation}
(here ${\bf Z}$ denotes the integer numbers), with coefficients
given by
\begin{equation}
 \varepsilon_n (\theta,\phi;\theta_0) = {1\over 2}\, a_n
  (\theta;\theta_0)\, \left(\cos n\phi +i\sin n\phi\right)\; .
 \label{14}
\end{equation}
The first two $a$ functions read
\begin{eqnarray}
 a_0(\theta;\theta_0) &=& {2\over \pi}\left(\phi_\star
  \cos\theta\cos\theta_0 + \sin\phi_\star \sin\theta\sin\theta_0
 \right)\; , \label{15}\\
 a_1(\theta;\theta_0)&=&{1\over \pi}\,\Bigl[2\sin\phi_\star\cos\theta
  \cos\theta_0 \nonumber \\ & &\quad +\left(\phi_\star+
  \sin\phi_\star\cos\phi_\star\right)\sin\theta\sin\theta_0\Bigr]
  \; . \label{16}
\end{eqnarray}
The auxiliary angle $\phi_\star$ is defined by
\begin{eqnarray}
 \cos\phi_\star(\theta,\theta_0) &=& -1  \qquad\qquad\qquad
  \mbox{\rm for }\theta< \theta_- \nonumber \\
 &=& -{\rm ctg}\theta\,{\rm ctg}\theta_0 \qquad\, \mbox{\rm for } \theta
  \in (\theta_-,\theta_+)  \label{17} \\
 &=& 1 \qquad\qquad\qquad\;\;\, \mbox{\rm for }\theta >\theta_+
  \nonumber
\end{eqnarray}
with $\theta_\pm = {\pi\over 2} \pm \theta_0$.
Similar expressions can also be found for higher--order terms,
but these will not be needed further on.
 
\subsubsection{Regular solution satisfying the boundary constraint}
After having discussed the source term, we are ready to find
a solution of the problem (\ref{9}) with the surface constraint
(\ref{10}). 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 a
convenient separation of the variables. We may write
\begin{equation}
 \Delta T'(x';\theta,\phi;\zeta) = \sum_{n\in{\bf Z}}\,
  t_{n}(x';\theta,\phi)\,\zeta^n \; . \label{18}
\end{equation}
Substituting this Fourier representation into (\ref{9}), we obtain a
system of decoupled second--order differential equations for the
amplitudes $t_n$. The general solution, which is regular inside the
body (i.e. for $x'\rightarrow \infty$), reads
\begin{equation}
 t_n = c_n\, \exp\left(-\sqrt{-in} x'\right)\; ,\label{19}
\end{equation}
with the constants $c_n$ to be determined by the boundary constraint
(\ref{10}). Easy algebra leads to
\begin{eqnarray}
 c_0 &=& {1\over\sqrt{2}}\, \left(\varepsilon_0 - \frac{1}{4}\right)
  \; \label{20} \\
 c_n &=& {\varepsilon_n\over \sqrt{2}}\, {(1+\lambda_n)-i\lambda_n
  \over 1+2\lambda_n+2\lambda_n^2} \equiv {\varepsilon_n\over \sqrt{2}}\,
  g_n \exp\left(i\delta_n\right)\; \label{21}
\end{eqnarray}
where $\lambda_n \equiv\Theta\sqrt{n}/2$.
 
\subsubsection{Thermal force and related quantities}
Having determined the surface temperature for each of the
elements we can now compute the thermal recoil force. We shall
assume Lambert's isotropic thermal emission geometry, as before, and
linearize
the fourth--power emission law. Calling ${\bf f}$ the thermal force
per unit of the body mass $m$ we obtain
\begin{equation}
 {\bf f} = -{2\sqrt{2}\over 3\pi}\,\alpha\Phi\,\int d\Omega\,\Delta T'
  \,{\bf n}\; , \label{22}
\end{equation}
where $\Phi = ({\cal E}_\star \pi R^2)/(mc)$ is the usual radiation
force factor.  The integration is performed on the whole sphere and
$d\Omega = \sin\theta d\theta d\phi$.  Finally, the resulting force is
expressed into a non--rotating reference frame $XYZ$, whose $Z$--axis
coincides with the body's spin axis and the $X$--axis coincides with
the $x$--axis of the rotating frame at time $t=0$.  The force
components in this new system read
\begin{eqnarray}
 f_X + if_Y &=& -{4\over 9}\,\alpha\,\Phi\,\sin\theta_0\,
  {(1+\lambda_1)+i\lambda_1\over 1+2\lambda_1+2\lambda_1^2}\; ,
  \label{23} \\
 f_Z &=& -{4\over 9}\, \alpha\, \Phi\, \cos\theta_0\; . \label{24}
\end{eqnarray}
This solution can be directly compared with Eqs.~(33,34) in
Vokrouhlick\'y (1998).  As expected, the two solutions coincide in
limit of ``large bodies'' ($R'\rightarrow \infty$ in Vokrouhlick\'y
1998).  In this approximation the spin--axis aligned component $f_Z$
is simply caused by diffusion of sunlight from the body, without any
influence of heat conduction.  On the contrary, the equatorial
components $f_X$ and $f_Y$ are influenced by heat conduction in the
body.  Our solution also matches well the solution of Peterson (1976).
However, we did not need to compute the involved Pad\'e representation
of the final force.  As mentioned earlier, the global choice of the
mean temperature $T_a$ on the sphere is the main source of this
simplification.
 
\subsection{Spheroidal bodies rotating around their axis of symmetry}
\subsubsection{Basics of spheroidal geometry}
Before we discuss the generalization of the previous solution to
the case of bodies with a more complex shape, we briefly recall some
basic results on the geometry of quadratic surfaces.
 
Consider an ellipsoid of rotation (hereinafter called a spheroid) with
equatorial radius $R$ and polar radius $(eR)$.  The parameter $e$
represents the geometric flattening of the body: $e\in (0,1)$ for
oblate spheroids and $e >1$ for prolate spheroids.  We need not to
introduce ellipsoidal coordinates.  Rather, we keep the
parametrization of the surface elements by the colatitude $\theta$
(measured from the axis of symmetry) and the longitude angle $\phi$
(measured in the equatorial plane from an arbitrary origin).  The
surface element $dS$ on the spheroid is given by
\begin{equation}
 dS = e^2 R^2 {J_4(\theta)\over J_2^4(\theta)}\,
  \sin\theta d\theta d\phi\; , \label{25}
\end{equation}
with the auxiliary functions of $\theta$ $(n=1,2,3, \ldots$):
\begin{equation}
 J_n(\theta) = \sqrt{e^n\sin^2\theta+\cos^2\theta}\; .
  \label{26}
\end{equation}
The total surface area of the spheroid is then given by $S = 4\pi
R^2\, s(e)$, where
\begin{equation}
 s(e) = \frac{1}{2}\left[1+{e^2 \over 2\eta}
  \ln\left({1+\eta\over 1-\eta}\right)\right]\; . \label{27}
\end{equation}
Here, we defined $\eta=\sqrt{1-e^2}$. Note that in the case of prolate
spheroids ($e>1$), $\eta$ becomes imaginary; in this case we can use
the following identity, holding for any complex number $z$:
\begin{equation}
 \frac{1}{2}\ln[(1+z)/(1-z)] = \frac{1}{i} {\rm arctg}(iz)\; .
  \label{28}
\end{equation}
A unit vector normal to the surface element at colatitude $\theta$
and longitude $\phi$ can be written as
\begin{equation}
  {\bf n}_\perp(\theta,\phi) = {1\over J_4(\theta)}\,
   \pmatrix{e^2\sin\theta\cos\phi \cr
            e^2\sin\theta\sin\phi \cr
                      \cos\theta \cr} \; . \label{29}
\end{equation}
Finally, we will also need the  geometric cross--section $P$ of the
spheroid when seen from a direction at an angle $\theta_0$
from the axis of symmetry. Simple algebra yields
\begin{equation}
 P(\theta_0) = \pi R^2 \, J_2(\theta_0)\; . \label{30}
\end{equation}
 
\subsubsection{Heat conduction problem in a spheroid}
Let us first describe the body's rotation.  Contrary to the uniform
rotation around a fixed axis which can be assumed for a spherical
body, the free rotation of an axisymmetric object is more complicated.
Generally, it consists of a proper rotation around the axis of
symmetry which, in turn, precesses around a constant direction in
space.  The thermal effects on a body whose axis of symmetry precesses
in space will be discussed in Sec.~2.3, and for the moment we shall
deal with the simpler case of a uniform rotation around the axis of
symmetry.  This will allow us to solve the problem analytically.  In
fact, most of the previous calculations remain unchanged.  We shall
just comment on a few differences.
 
First, we have the definition of the mean temperature $T_a$. As
before, this quantity follows from the balance of the absorbed and
re--radiated energy over one rotation cycle of the body. Using the
variables defined above we have
\begin{equation}
 -\alpha{\cal E}_\star\, \int_{|\zeta|=1}{id\zeta\over \zeta}P(\theta_0)
  = S \epsilon \sigma T^4_a\; . \label{31}
\end{equation}
Scaling $T_a$ by the reference temperature $T_\star$ we obtain
\begin{equation}
 T'_a = {1\over \sqrt{2}}\,\left[{J_2\left(\theta_0\right) \over
  s\left(e\right)}\right]^{1/4}\equiv{1\over \sqrt{2}}\, \mu(\theta_0;e)
  \; . \label{32}
\end{equation}
Given this expression for the mean temperature, we note that the
linearized boundary constraint (\ref{10}) for the temperature
variation $\Delta T'$ reads
\begin{equation}
 \sqrt{2}\mu^3\Delta T'- \Theta\left[{\partial \Delta T'\over \partial
  x'}\right]_0 = \alpha\,\Delta {\cal E}'\; ,\label{33}
\end{equation}
with the source term $\Delta {\cal E}'={\cal E}'-\frac{1}{4}\mu^4$.
The Fourier representation (\ref{13}) of this source term is
unchanged, but the coefficients are now given by
\begin{equation}
 \varepsilon_n (\theta,\phi;\theta_0) = {1\over 2}\, {a_n
  (\theta;\theta_0)\over J_4(\theta)}\, \left(\cos n\phi +i\sin n\phi
  \right)\; . \label{34}
\end{equation}
The first two functions $a(\theta;\theta_0)$ now read
\begin{eqnarray}
 a_0(\theta;\theta_0) &=& {2\over \pi}\left(\phi_\star
  \cos\theta\cos\theta_0 + e^2\sin\phi_\star \sin\theta\sin\theta_0
 \right)\; , \label{15bis}\\
 a_1(\theta;\theta_0)&=&{1\over \pi}\,\Bigl[2\sin\phi_\star\cos\theta
  \cos\theta_0 \nonumber \\ & &\quad +e^2\left(\phi_\star+
  \sin\phi_\star\cos\phi_\star\right)\sin\theta\sin\theta_0\Bigr]
  \; . \label{16bis}
\end{eqnarray}
and the auxiliary angle $\phi_\star$ is defined by
\begin{eqnarray}
 \cos\phi_\star(\theta,\theta_0) &=& -1  \qquad\qquad\qquad\quad
  \mbox{\rm for }\theta< \theta_- \nonumber \\
 &=& -e^{-2}\,{\rm ctg}\theta\,{\rm ctg}\theta_0 \qquad \mbox{\rm for }
  \theta\in (\theta_-,\theta_+)  \label{35} \\
 &=& 1 \qquad\qquad\qquad\qquad\; \mbox{\rm for }\theta >\theta_+
  \nonumber
\end{eqnarray}
with ${\rm ctg}\theta_\pm = \mp e^2{\rm tg}\theta_0$.
 
Because the heat conduction equation (\ref{9}) is independent of the
particular geometry of the body, the Fourier expansion (\ref{18}) with
coefficients given by (\ref{19}) is still valid. The only difference
is that the integration constants $c_n$ now must satisfy
(\ref{33}), which yields
\begin{eqnarray}
 c_0 &=& {1\over\sqrt{2}\mu^3}\, \left(\varepsilon_0 - \frac{1}{4}\mu^4
  \right) \label{36} \\
 c_n &=& {\varepsilon_n\over \sqrt{2}\mu^3}\, {(1+\lambda_n)-i\lambda_n
  \over 1+2\lambda_n+2\lambda_n^2} \equiv {\varepsilon_n\over \sqrt{2}\mu^3}
  \, g_n \exp\left(i\delta_n\right)\; , \label{37}
\end{eqnarray}
with $\lambda_n \equiv\Theta\sqrt{n}/(2\mu^3)$.
 
\subsubsection{Diurnal Yarkovsky force on a spheroid}
Given the solution for the surface temperature derived above,
we can now compute the total thermal recoil acceleration. Integrating
over the surface of the spheroid (with $d\Omega = \sin\theta\, d\theta
d\phi$ as before), we have
\begin{equation}
 {\bf f} = -{2\sqrt{2}\over 3\pi}\,\alpha\Phi\,e^2\mu^3\,\int
  d\Omega\,{J_4(\theta)\over J_2^4(\theta)}\,\Delta T'\,{\bf
  n}_\perp\; . \label{38}
\end{equation}
Interestingly, this two--dimensional integration can be worked out
analytically
as in the case of a spherical body. The thermal acceleration,
projected onto the unit vectors of the non--rotating $XYZ$ system,
has components:
\begin{eqnarray}
 f_X + if_Y &=& -{4\over 9}\,\alpha\,\Phi\,\sin\theta_0\,
  {(1+\lambda_1)+i\lambda_1\over 1+2\lambda_1+2\lambda_1^2}\,
  \psi_{XY}(e)\; , \label{39} \\
 f_Z &=& -{4\over 9}\, \alpha\, \Phi\, \cos\theta_0\, \psi_Z(e)
  \; . \label{40}
\end{eqnarray}
We recall that
\begin{equation}
 \lambda_1 = {\Theta\over 2}\, \left[{s\left(e\right)
  \over J_2\left(\theta_0\right)}\right]^{3/4}\; , \label{41}
\end{equation}
and the two functions of $e$ in the right--hand sides
of (\ref{39}) and (\ref{40}) read
\begin{eqnarray}
 \psi_{XY}(e) &=& {3\over 4} {e^2\over \eta^2}\left[{1+\eta^2\over
  2\eta}\ln\left({1+\eta\over 1-\eta}\right)-1\right]\; , \label{42}\\
 \psi_Z(e) &=& {3\over 2\eta^2}\left[1-{e^2\over 2\eta}
  \ln\left({1+\eta\over 1-\eta}\right)\right] \; . \label{43}
\end{eqnarray}
 
The two functions $\psi$ have been plotted in Figure~\ref{f1}.  The
equatorial $X$ and $Y$ force components are reduced by the
corresponding $\psi_{XY}$ factors in the case of oblate spheroids and
{\it vice versa}.  This result is consistent with the intuitive idea
that the normal unit vectors ${\bf n}_\perp$ to the polar surface
elements of oblate spheroids are more aligned with the $Z$--axis than
the corresponding normal vectors on a sphere.  On the contrary, the
${\bf n}_\perp$ vectors to the equatorial surface elements lean toward
the equatorial plane for prolate spheroids, and therefore contribute
to enhancing the equatorial components of the thermal force.
 
\begin{figure}[t]
\epsfxsize=8.7cm
\begin{center}
\mbox{\epsfbox{aa14f1.eps}} \end{center} \vspace{-2.2cm}
\caption[]{Non--sphericity factors $\psi_{XY}$ and $\psi_Z$ vs. the
flattening parameter $e$.
\label{f1}}
\end{figure}
 
However, when we try to decide whether the thermal force
components (\ref{39}) and (\ref{40}) are smaller or greater than
those on a sphere, we should take into account two more factors:
\begin{itemize}
\item the appearance of the $\mu(\theta_0;e)$ factor in the
      $\lambda_1$ parameter;
\item the equatorial radius $R$ of the spheroid appearing in
      the definition of the $\Phi$ radiation force factor.
\end{itemize}
The latter item is related to the problem of the choice of the sphere
whose thermal effects should be compared with those for the spheroid.
Of course, different such choices are possible.  Hereinafter, bodies
having the same mass will be compared.  Provided the mean density is
factorized, equal masses imply that the radius of the comparison
sphere has to be $(e^{1/3}\, R)$. Supposing that the thermal parameter
$\Theta$ is sufficiently large, we can approximate the rational
function of $\lambda_1$ in (\ref{39}) by $\frac{1}{2}(1+i)/
\lambda_1$.  Then, we conclude that the equatorial thermal force
components $f_X$ and $f_Y$ from (\ref{39}) can be written as the
corresponding force components on a sphere (\ref{23}) multiplied by a
numerical factor $\xi_{XY}$.  In a similar way, we define a factor
$\xi_Z$ for the third force component $f_Z$.  Putting together the
previous results, we have
\begin{eqnarray}
 \xi_{XY} &\simeq& e^{-2/3}\,\psi_{XY}(e)\,
  \left[J_2\left(\theta_0\right)\over {s\left(e\right)}
  \right]^{3/4}\; , \label{44}\\
 \xi_Z &=& e^{-2/3}\,\psi_Z(e)\; . \label{45}
\end{eqnarray}
\begin{figure}[t]
\epsfxsize=8.7cm
\begin{center}
\mbox{\epsfbox{aa14f2.eps}}
\end{center}
\vspace{-0.5cm}
\caption[]{Curves corresponding to a constant $\xi_{XY}$ factor
projected on the plane defined by the solar tilt angle $\theta_0$
(smaller than $45^\circ$) and the body's flattening parameter $e$.
The values of $\xi_{XY}$ label the different curves: $\xi_{XY}>1$
means that the thermal equatorial components are larger than for a
spherical body with the same mass and {\it vice versa}.
\label{f2}}
\end{figure}
 
\begin{figure}[t]
\epsfxsize=8.7cm
\begin{center}
\mbox{\epsfbox{aa14f3.eps}}
\end{center}
\vspace{-0.5cm}
\caption[]{The same as in Figure 2 but for $\theta_0$ greater
 than $45^\circ$.
 \label{f3}}
\end{figure}
 
Figures~\ref{f2} and \ref{f3} show contour plots of the $\xi_{XY}$
factor in the parameter plane $\theta_0$ vs. the flattening parameter
$e$.  We note that the maximum differences in the equatorial force
components occur for $\theta_0 = 90^\circ$, when the diurnal effect
itself is maximum [because of $\sin\theta_0$ factor in
Eq.~(\ref{39})].  The greatest mismodelling (by a factor $2$ to $3$)
is observed when $e$ reaches values as low as $0.3$ or as high as $3$.
The $\xi_Z$ factor has the same dependence on $e$ as the previously
discussed function $\psi_Z$ --- see Figure~\ref{f4}.
 
\begin{figure}[t]
\epsfxsize=8.7cm
\begin{center}
\mbox{\epsfbox{aa14f4.eps}}
\end{center}
\vspace{-2.2cm}
\caption[]{$\xi_Z$ vs. the flattening parameter $e$.
 \label{f4}}
\end{figure}
 
\subsection{Spheroidal bodies rotating around an arbitrary axis}
In this section we shall generalize the previous results to the case
when the symmetry axis of the spheroid tumbles around a fixed
direction in space.  Such a free--precessional state seems to be a
natural rotational mode for meter-sized asteroidal fragments, since
the internal dissipation of energy is not efficient enough to align
the spin axis with the angular momentum vector in a short time.
According to the theory of Burns \& Safronov (1973) the typical
time scale for such a dissipation process is about $2$~Byr
(considering rather conservative values of the dissipation parameter),
whereas in the main asteroid belt collisions transferring relatively
large amounts of angular momentum typically occur every few Myr only,
and for stony bodies shattering impacts occur every few tens of Myr
(Farinella et~al. 1998).  A generalization of the previous results to
the tumbling rotation case appears also important for another reason:
as the spin axis sweeps the precession cone, different parts of the
surface are more efficiently illuminated by the radiation flux.  One
might expect that a more uniform distribution of surface temperatures
will be established, decreasing the amplitude of the thermal
acceleration.
 
The orientation in space of the reference frame rigidly rotating with
the body (whose $z$-axis coincides with the symmetry axis) will be
described by three Euler angles:
\begin{itemize}
\item $\psi_E = \omega t$ is the precession angle describing the
 motion of the  projection of the symmetry axis on the $XY$ equatorial
 plane of the non--rotating reference frame, whose
 $Z$-axis coincides with the angular momentum vector;
\item $\theta_E = \tau=$~constant is the nutation angle between the
 symmetry axis ($z$) and the $Z$-axis of the non--rotating system;
\item $\varphi_E$ is the angle of proper rotation.
\end{itemize}
To be consistent with the previous notations, we shall describe the
transformation between the two frames by a 3--2--3 sequence of
individual rotations, meaning that the second rotation (by the
constant nutation angle $\theta_E$) is performed around the $y$-axis.
The Euler angles used in this paper thus differ slightly from the
classical set of Euler angles given in the 3--1--3 sequence.
 
Since the precession frequency is $\omega$, taking $\zeta=
\exp(i\omega t)$ as before, the proper rotation velocity is
\begin{equation}
 {1\over \omega}{d\varphi_E\over dt} = -{C-A\over C}\,\cos\tau =
  -\frac{\eta^2}{2}\,\cos\tau \label{46}
\end{equation}
(see e.g.  Landau \& Lifschitz 1960), where $C$ and $A$ are the
principal moments of inertia around the axis of symmetry and an
equatorial axis, respectively, and for homogeneous spheroids we have
$(C-A)/C = \frac{1}{2} \,\eta^2$.
 
The general scheme outlined in Sec.~2.1 is unchanged, because it
concerns the individual thermal histories of different surface
elements of the body.  What becomes more complicated is mainly the
definition of the mean temperature $T'_a$ and the radiation flux term
$\Delta{\cal E}'$.  We shall now discuss these quantities in some
detail.
 
Consider a surface element $dS$ specified by spherical angles
$\theta$ and $\phi$ in the system rigidly rotating with the body.
The external radiation flux through this element (scaled by
${\cal E}_\star$) is
\begin{eqnarray}
 {\cal E}'(\theta,\phi;\theta_0,\zeta) &=& {\bf R}^2_\tau\,
  {\bf R}^3_{\varphi_E(\zeta)}\,{\bf n}_\perp(\theta,\phi)
  \cdot{\bf n}_0(\theta_0,\zeta)\nonumber \\
  & & \qquad\qquad\qquad\qquad\quad
   \mbox{\rm if not negative,}  \nonumber \\
  &=& 0 \qquad\qquad\qquad\qquad\;\; \mbox{\rm otherwise,}
  \label{47}
\end{eqnarray}
 
which generalizes Eq.~(\ref{12}), valid for a spherical body. ${\bf
R}^2_\tau$ and ${\bf R}^3_{\varphi_E}$ are the usual orthonormal
rotation matrixes around the $y-$ and $z$-axes (the angle is specified
by the index).  The simplicity of the solution for thermal effects
from the previous Secs.~2.1 and 2.2 is mainly based on the linearity
of the heat conduction equation (\ref{9}) and the boundary constraints
(\ref{10}) or (\ref{33}).  However, the exact periodicity of the
solution, expressed by the Fourier expansions of all the important
variables [in particular the temperature $T'$ in Eq.~(\ref{18})] was
an equally important assumption.  The ortonormality of the Fourier
modes allowed us to decouple the differential equations for the $t_n$
amplitudes of the temperature terms in Eq.~(\ref{18}).  This second
key property of the diurnal thermal effect solution is violated in the
case of a tumbling spin axis, because the precession frequency ${\dot
\psi}_E=\omega$ and the proper rotation frequency ${\dot \varphi}_E=
-\frac{1}{2} \eta^2 \cos\tau\, \omega$ -- in general are not
commensurable with each other.
 
In order to keep our solution simple enough, we shall adopt a
two--steps simplifying approach:
\begin{itemize}
\item first, we shall neglect the proper rotation of the body,
 assuming $\frac{1}{2} \eta^2 \cos\tau \ll 1$; then the rotation
 matrix $R^3_{\varphi_E(\zeta)}$ can be approximated by a unit,
 time independent matrix;
\item second, we shall verify the validity of the
 above--mentioned solution on resonance ``slices'' defined by ${\dot
 \varphi}_E/ {\dot \psi}_E = -\frac{1}{2} \eta^2 \cos\tau = -1/n$,
where $n$ is an integer number $\geq 1$.
\end{itemize}
 
In the latter case the periodicity of the diurnal effect solution is
that of the proper rotation ($n$ times the precession period).
Although we shall not solve the problem in the most general case, we
believe that our results are rich enough to give an idea about the
role of non--sphericity effects for tumbling fragments.
 
\subsubsection{Proper rotation neglected}
In the first step of our approximation we disregard the
proper rotation of the body around its symmetry axis. We can then
expand the insolation function $\Delta {\cal E}'$ of the boundary
condition (\ref{33}) in the Fourier series (\ref{13}), with
coefficients $\varepsilon_n$ computed numerically (for $n\neq 0$) by
\begin{equation}
 \varepsilon_n(\theta,\phi;\theta_0) = {1\over 2\pi i}\, \int_{|\zeta|=1}
  \,{d\zeta\over \zeta^{n+1}}\,{\cal E}'(\theta,\phi;\theta_0,\zeta)
  \; . \label{48}
\end{equation}
The $(n=0)$ coefficient constrains, like in the previous
discussion, the value of the mean temperature $T'_a\equiv\mu(\theta_0;
e,\tau)/ \sqrt{2}$. The balance between the total energy absorbed
by the body over one precessional cycle and the amount of thermally
re--radiated  energy over the same period is expressed by
\begin{equation}
 \int d(\cos\theta)d\phi\, {J_4(\theta)\over J_2^4(\theta)}\,
  \varepsilon_0(\theta,\phi;\theta_0)=0\; , \label{49}
\end{equation}
with the integration performed over the whole surface of the body and
$\varepsilon_0(\theta,\phi;\theta_0)$ given by
\begin{eqnarray}
 \varepsilon_0(\theta,\phi;\theta_0) &=& {1\over 2\pi i}\,
  \int_{|\zeta|=1}\,{d\zeta\over \zeta}\,\left[{\cal E}'
  (\theta,\phi;\theta_0,\zeta)-\frac{1}{4}\mu^4\right]\; .
  \nonumber \\& & \label{50}
\end{eqnarray}
Eqs.~(\ref{49}) and (\ref{50}) provide the $\mu$ factor appearing in
the mean temperature in the form
\begin{equation}
 \mu\left(\theta_0;e,\tau\right) = \left[{{\cal I}(\tau,\theta_0;e)
  \over s(e)}\right]^{1/4} \; , \label{51}
\end{equation}
where the integral ${\cal I}$ is given by
\begin{equation}
 {\cal I}(\tau,\theta_0;e)={e\over \pi}
  \int_{c_+}^{c_-}dx\,\sqrt{{1+\nu x^2\over (x-c_+)(c_- -x)}}\; .
 \label{52}
\end{equation}
Here, we denoted $\nu=(1-e^2)/ e^2$ and $c_\pm = \cos(\theta_0\pm
\tau)$. In general, the right--hand side of (\ref{52}) should be
computed numerically, but in some particular cases it is integrable
analytically. For later use, we note that in the case of the
configuration $\theta_0=90^\circ$ (solar direction normal to the
rotation axis) we obtain
\begin{eqnarray}
 {\cal I}(\tau,\pi/2;e)&=&  {2e\over \pi}\,
  \sqrt{1+\nu\sin^2\tau}\, E\left(\sqrt{{\nu\sin^2\tau\over 1+\nu
  \sin^2\tau}}\right)\nonumber \\& & \qquad\qquad\qquad\qquad\quad
   \quad {\rm for }\;\;e\leq 1\; ,\label{53}\\
  &=& {2e\over \pi}\, E\left(\sqrt{-\nu\sin^2\tau}\right)
  \quad\;\, {\rm for }\;\;e\geq 1\; .
\end{eqnarray}
Here $E$ is the complete elliptic integral of the second kind.
Expressions for $e\leq 1$ and $e\geq 1$ can be obtained from each
other by means of the complex argument transformation of elliptic
functions:
\begin{equation}
 E\left(i{k\over\sqrt{1-k^2}}\right) = {E(k)\over\sqrt{1-k^2}}
  \label{53bis}
\end{equation}
(see e.g.  Byrd \& Friedman 1971).  Having pre--computed numerically
all the previous functions, we can determine the thermal force
averaged over one precessional cycle in the non--rotating inertial
system.  The temperature of each surface element can be expanded in
the Fourier series (\ref{18}), with the amplitude functions given by
(\ref{19}).  The surface boundary constraint then yields the
integration constants (\ref{36}) and (\ref{37}).  Finally, the thermal
force is obtained by numerical integration of (\ref{38}) over the
whole surface of the spheroid, with the unit normal vector ${\bf
n}_\perp$ of (\ref{29}) transformed into the non-rotating system
using the rotation matrix ${\bf R}^3_{\psi_E}\,{\bf R}^2_\tau$.  After
performing all the necessary algebra we obtain
\begin{eqnarray}
 f_X &=& -{2\over 3}\alpha\, \Phi\,{e^2\over 1+2\lambda_1+
  2\lambda_1^2}\times\nonumber\\
  & &\qquad\quad \Bigl\{e^2\cos\tau\left[\left(1+\lambda_1
  \right){\rm Re}E_{1,X}+\lambda_1{\rm Im}E_{1,X}\right]\nonumber\\
 & & \qquad\quad +\sin\tau\left[\left(1+\lambda_1\right)
  {\rm Re}E_{1,Z}+\lambda_1{\rm Im}E_{1,Z}\right]\nonumber\\
 & & \qquad\quad -e^2\left[\lambda_1{\rm Re}E_{1,Y}-
  \left(1+\lambda_1\right){\rm Im}E_{1,Y}\right]\Bigr\}\; ,\label{54}\\
 f_Y &=& -{2\over 3}\alpha\, \Phi\,{e^2\over 1+2\lambda_1+
  2\lambda_1^2}\times\nonumber\\ & &
  \qquad\quad  \Bigl\{e^2\cos\tau\left[\lambda_1{\rm Re}
  E_{1,X}-\left(1+\lambda_1\right){\rm Im}E_{1,X}\right]\nonumber\\
 & & \qquad\quad +\sin\tau\left[\lambda_1{\rm Re}E_{1,Z}-
  \left(1+\lambda_1\right){\rm Im}E_{1,Z}\right]\nonumber\\
 & & \qquad\quad +e^2\left[\left(1+\lambda_1\right)
  {\rm Re}E_{1,Y}+\lambda_1{\rm Im}E_{1,Y}\right]\Bigr\}\; ,\label{55}\\
 f_Z &=& -{2\over 3}\alpha\, \Phi\,e^2
  \left(\cos\tau E_{0,Z}-e^2\sin\tau E_{0,X}\right)
  \label{56}
\end{eqnarray}
for the three thermal force components. The auxiliary vector
quantities ${\bf E}_n$ are defined by
\begin{equation}
 {\bf E}_n(\theta_0) = \int d(\cos\theta)d\phi\, {J_4(\theta)
  \over J_2^4(\theta)}\,\varepsilon_n(\theta,\phi;\theta_0)\,
  {\bf n}_\perp\; ,  \label{57}
\end{equation}
and the thermal parameter coefficient $\lambda_1$ reads
\begin{equation}
 \lambda_1 = {\Theta\over 2}\, \left[{s(e)\over
  {\cal I}(\tau,\theta_0;e)}\right]^{3/4}\; . \label{58}
\end{equation}
Obviously, Eqs.~(\ref{54}) -- (\ref{56}) become equal to (\ref{39})
and (\ref{40}) for $\tau=0$ (so that the symmetry axis of the body
coincides with the rotation axis).  This property can be easily
verified since in this limit ${\rm Re}(E_{1,X}) = {\rm Im}(E_{1,Y})$
and ${\rm Re}(E_{1,Y}) = {\rm Im}(E_{1,X})$.
\begin{figure}[t]
\epsfxsize=8.7cm
\begin{center}
\mbox{\epsfbox{aa14f5.eps}}
\end{center}
\vspace{-0.5cm}
\caption[]{Contour plot of the $\xi_{XY}$ function in the plane
 defined by the spin axis tilt $\tau$ from the axis of symmetry and
 the flattening parameter. Each contour line is labelled by the
 corresponding value of $\xi_{XY}$. The solar direction is
 assumed to be normal to the spin axis
 ($\theta_0=90^\circ$). In the shaded area the proper rotation
 frequency ${\dot \varphi}_E$ is greater than $\frac{1}{10}\,{\dot
 \psi}_E$.
 \label{f5}}
\end{figure}

Adopting the methodology of Sec.~2.1 to compare the thermal force
components (\ref{54}) -- (\ref{56}) with the corresponding thermal
force acting on a sphere, we define the ``correction factors''
$\xi_{XY}$ and $\xi_Z$ by: $\xi_{XY} = f_X/f^{\rm sphere}_X =
f_Y/f^{\rm sphere}_Y$ and $\xi_Z =f_Z/f^{\rm sphere}_Z$.  A
sufficiently large value of the thermal parameter $\Theta$ is assumed
and bodies having the same total mass are compared.  Figure~\ref{f5}
shows contour curves of $\xi_{XY}$ in the plane $\tau$ vs. $e$.  As an
example, we have assumed the special configuration $\theta_0=90^\circ$
in which the thermal equatorial force components are maximum.  The
results for $\tau=0$ match those for $\theta_0= 90^\circ$ in
Figure~\ref{f3}.  A reversal in the magnitude of the equatorial force
components compared to those acting on a sphere at large values of
$\tau$ has to be noted.  For small values of $\tau$ and oblate
spheroids ($e<1$), for instance, the surface normal vectors ${\bf
n}_\perp$ in the integral (\ref{38}) are preferentially directed along
the $Z$-axis, thus decreasing the amplitude of the equatorial force
components.  The corresponding value of the $\xi_{XY}$ factor is
smaller than unity.  The opposite occurs for values of $\tau \simeq
90^\circ$ (symmetry axis perpendicular to the rotation axis): in this
case the surface normal vectors ${\bf n}_\perp$ are preferentially
tilted toward the equatorial plane, thus contributing efficiently to
the $(f_X,f_Y)$ force components.  The corresponding $\xi_{XY}$ factor
is larger than unity.  Obviously, the situation is just opposite in
the case of prolate spheroids ($e>1$).
 
The shaded areas in Figure~\ref{f5} correspond to the configurations
where the ratio between the proper rotation frequency and the
precession frequency ${\dot \varphi}_E/{\dot \psi}_E =
-\frac{1}{2}\eta^2\, \cos\tau$ is larger in absolute magnitude than
$\frac{1}{10}$.  In this region our simplified solution, neglecting
the proper rotation of the body, may not be justified.  The large
value of the thermal parameter $\Theta$ means that the thermal memory
time scale is much longer than the rotation period.  In this
situation, one might expect that the proper rotation could contribute
to smearing the surface temperature, decreasing the amplitude of the
thermal effects.  In the next paragraphs we assess briefly the
importance of such phenomena.
 
\subsubsection{The influence of the proper rotation}
In the most general case, when the proper rotation angle $\varphi_E$
is variable, we note that the scaled radiation flux
(\ref{49}) may not be a periodic function. This would violate
the basic assumptions of our solution, as discussed above.
In order to keep ${\cal E}'$ periodic, we shall investigate the
thermal effects in the ``resonant'' cases: ${\dot \varphi}_E/{\dot
\psi}_E = -1/n$, where $n$ is an integer. We believe that
these cases are representative enough to draw general conclusions.
 
In the resonant situations the periodicity of ${\cal E}'$ and
of the whole solution is recovered. The basic time scale is now given
by the proper rotation cycle, over which the body undergoes $n$
precessional cycles. All the necessary formul{\ae} for computing the
thermal force have been given above. However, as the algebra is rather
lengthy, we shall not give here the final results for the
thermal force components explicitly.
\begin{figure}[t]
\epsfxsize=8.7cm
\begin{center}
\mbox{\epsfbox{aa14f6.eps}}
\end{center}
\vspace{-2.2cm}
\caption[]{$\xi_{XY}^n$ vs. the flattening parameter $e$ for some
$n:1$ ``resonance slices''.  The degree $n$ of the resonance labels
the curves.  The dashed lines correspond to the simple approximation
when the proper rotation of the body around its symmetry axis is
neglected (results from Figure~5), while the solid lines account for
the proper rotation.
\label{f6}}
\end{figure}
 
Figure~\ref{f6} shows the $\xi_{XY}^n$ factor vs. the flattening
parameter $e$ for three resonance conditions: $n=2,5,10$.  The dashed
lines are replotted from Figure~\ref{f5} to allow a comparison with
the case when the proper rotation is neglected, while the solid lines
correspond to the complete solution.  At the $1:10$ resonance there is
a very good agreement between the two results.  At the next ``slice'',
corresponding to the $1:5$ resonance, the agreement is still
acceptable, whereas in the case of the $1:2$ resonance (precession
frequency twice as large as the proper rotation frequency) we observe
a significant difference.  Interestingly, in the exact solution the
force components are larger than the corresponding components from the
simplified solution.
 
\begin{figure}[t]
\epsfxsize=8.7cm
\begin{center}
\mbox{\epsfbox{aa14f7.eps}}
\end{center}
\vspace{-2.2cm}
\caption[]{Mean value of the $\xi_{XY}$-factor over a random sampling
of the nutation angle $\tau$ and the solar tilt angle $\theta_0$
($500$ values considered), as a function of the flattening parameter
$e$.  The proper rotation of the body is neglected.  Error bars
correspond to the Gaussian standard deviation.
 \label{f7}}
\end{figure}
 
\subsubsection{Results averaged over the orientation of the spin
axis in space}
 
In the previous discussion we have compared the exact solution for the
diurnal Yarkovsky force on spheroids with the corresponding solution
for a sphere using a special configuration, that is to say, we had
fixed a particular value of the $\theta_0$ angle between the solar
direction and the body's spin axis and/or of the $\tau$ angle between
the symmetry and the spin axes.  In reality, these geometric
parameters change owing to the body's revolution around the radiation
source (the Sun) and the collisions with other asteroid fragments.  In
both cases the corresponding time scale is typically much shorter than
the period over which the orbital evolution is studied by numerical
integrations.  Therefore, we might be interested in averaging the
$\xi_{XY}$ factor over all the possible configurations of the
$(\theta_0$ and $\tau)$ angles.  The mean values of the $\xi_{XY}$
factor over a sample of $500$ random values of the angles together
with the corresponding Gaussian standard deviations are shown in
Figure~\ref{f7}.  In contrast to the previous results, we get a much
better agreement between the diurnal Yarkovsky force computed for
spheroids and for a sphere.  The peak differences are of the order of
$30\%$, but in most cases are smaller than $20\%$.  As expected, the
dispersion of the results --- indicated by the standard deviations ---
is larger for more elongated bodies.
 
\section{Conclusions}
The main results of this paper may be summarized as follows:
\begin{itemize}
 
\item The diurnal Yarkovsky recoil force has been computed for bodies
whose shape lacks a spherical symmetry, that is for axisimmetrical
spheroids with an arbitrary flattening.  The only important
simplifying assumption of our approach is that the average size of the
body must be much larger than the penetration depth of the thermal
diurnal wave.
\item The results are derived in a closed analytical form if the body
rotates around its axis of symmetry.  In the more general case of a
free--precessing object, the results are given numerically.
\item Our main goal has been that of comparing the diurnal Yarkovsky
force acting on an irregularly shaped body with that acting on a
sphere, since this shape corresponds to a commonly used approximation
in the long--term numerical integration schemes for asteroidal
fragment orbits.  Our results can be used to assess the corresponding
mismodelling of the Yarkovsky perturbations found from these
integrations.
\item Our results indicate that, depending on the geometric parameters
(degree of flattening, nutation angle etc.), the simple
spherical--shape approximation for the diurnal Yarkovsky force
provides results mismodelled up to a factor of $2$ or $3$.
\item However, when an average over all possible orientations of the
spin axis and the solar direction is performed, there is a much better
agreement between the thermal force acting on a sphere and our results
for spheroids.  This suggests that neglecting the influence of
non--sphericity effects may not degrade significantly the results of
long--term orbital integrations including the diurnal Yarkovsky
effect.
\end{itemize}
 
Of course, there is the possibility that the resulting thermal force
may be significantly different for bodies of much more irregular shape
than those studied in this paper.  This issue is beyond the scopes of
the current study and may be require further work in the future.
However, it seems very likely that the results will not differ too
much in the most general case of a tumbling spin axis, as discussed
above.
 
The analysis of this paper suggests that non--sphericity effects
should not be seen as a very serious concern in modelling the diurnal
Yarkovsky perturbations.  Therefore, in a forthcoming paper we plan to
apply the thermal model for the diurnal Yarkovsky effect on a sphere
as derived in paper~I to a long--term integration of fragment orbits.
Several important issues related to the delivery of meteorites and
asteroidal fragments to the inner Solar System will be considered.
 
\medskip
\noindent{\bf Acknowledgements}\hskip 2mm
The author thanks M.~Bro\v z for his help in drawing Fig.~5, W.F.~Bottke
and P.~Farinella for helpful discussions.  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 Binzel R.P., P. Farinella, V. Zappal\'a, A. Cellino, 1989, in:
 Asteroids~II, eds. R.P. Binzel, T. Gehrels and M.S. Matthews,
 Univ. of Arizona Press, Tucson, p.~416
\item Bottke W.F., D.P. Rubincam, J.A. Burns, Consequences of the
 Yarkovsky effect on the orbital evolution of meteoroids: A
 prospectus, paper presented at Lunar Planet. Sci. Conf., Houston,
 1998
\item Burns J.A., V.S. Safronov, 1973, MNRAS 165, 403
\item Byrd P.F., M.D. Friedman, 1971, Handbook of Elliptic Integrals
 for Engineers and Scientists, Springer--Verlag, Berlin
\item Catullo V., V. Zappal\'a, P. Farinella, P. Paolicchi, 1984,
 AA 138, 464
\item Farinella P., D. Vokrouhlick\'y, W.K. Hartmann, 1998, Icarus,
 in press
\item Fujiwara A., G. Kamimoto, A. Tsukamoto, 1978, Nature 272, 602
\item Fujiwara A., P. Cerroni, D.R. Davis, et~al., 1989,
 in: Asteroids~II, eds. R.P. Binzel, T. Gehrels and M.S. Matthews,
 Univ. Arizona of Press, Tucson, p.~240
\item Giblin I., G. Martelli, P.N. Smith et~al., 1994, Icarus 110,
 203
\item Giblin I., P. Farinella, 1997, Icarus 127, 424
\item Harris A.W., 1994, Icarus 107, 209
\item Hartmann W.K., P. Farinella, D. Vokrouhlick\'y, S.J.
 Weidenschilling, A. Morbidelli, F. Marzari, D.R. Davis, E. Ryan, 1998,
 Meteoritics Planet. Science, in press
\item Landau L.D., E.M. Lifschitz, 1960, Mechanics, Addison--Wesley,
 Reading
\item Ostro S.J., K.D. Rosema, R.S. Hudson, et~al., 1995,
 Nature 375, 474
\item Peterson C., 1976, Icarus 29, 91
\item Radzievskii V.V., 1952, Astron. Zh. 29, 162 (in Russian)
\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 Vokrouhlick\'y D., 1998, AA, in press (paper~I)
\item Vokrouhlick\'y D., P. Farinella, 1998, Astron. J., submitted
\item Zappal\'a V., P. Farinella, 1997, Adv. Space Res. 19, 181
\end{description}
\end{document}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
