%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% A two-layer model of the Yarkovsky seasonal effect of sphere
% ver. DV March 25 1999, corrected by PF Apr 1 1999
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\documentstyle[epsf]{l-aa}
%\renewcommand{\baselinestretch}{2}
%\onecolumn
\thesaurus{05.03.1; 07.13.2}
\title{An improved model of the seasonal Yarkovsky force
 for regolith--covered asteroid fragments}
\author{David Vokrouhlick\'y and Miroslav Bro\v z}
\institute{Institute of Astronomy, Charles University,
    V Hole\v{s}ovi\v{c}k\'ach 2, CZ-180 00 Prague 8, Czech Republic\\
    (Electronic addresses: vokrouhl@mbox.cesnet.cz,
    miroslav.broz@usa.net)}
\offprints{D. Vokrouhlick\'y}
\date{\today}
\begin{document}
\maketitle
\markboth{D. Vokrouhlick\'y and M. Bro\v z: An improved model of the
 Yarkovsky seasonal effect...}{}
\begin{abstract}
 We consider an exact solution of the seasonal Yarkovsky effect,
 the mean--motion frequency mode of the recoil force due to
 reradiated sunlight, on a spherical asteroid fragment. The body is
 assumed to have a low--conductivity (``regolith--like'') surface layer
 covering a more conductive core. The penetration depth of the
 seasonal thermal wave in the low--conductivity surface material is
 typically larger than its geometrical thickness, a fact that has
 not been taken into account consistently by the previous simplified
 solutions. Our approach removes this drawback and the results
 indicate (i) a significantly increased seasonal mobility of the
 dekameter--sized (and larger) fragments with an insulating surface
 layer and (ii) a significantly decreased seasonal mobility of the
 meter--sized fragments with this structure. These results may affect
 the accuracy of simulations of the  meteorite and NEA delivery to the
 Earth.
\keywords
{Celestial mechanics -- Minor planets, asteroids -- Meteoroids}
\end{abstract}

\section{Introduction}
Thermal phenomena on asteroid fragments have attracted considerable
attention over the last few years since their dynamical aspects, known as
the Yarkovsky effects, were found to be important in solving some
interesting problems in meteoritics and the delivery of NEAs to the Earth
vicinity (e.g. Farinella et~al. 1998, Hartmann et~al. 1999,
Vokrouhlick\'y \& Farinella 1998, Farinella \& Vokrouhlick\'y 1999).
A detailed understanding of the real astronomical role of the Yarkovsky
effects requires (i) a precise knowledge of the thermal and other
parameters affecting their efficiency (e.g. surface thermal conductivity
of the asteroid fragments, details of their collisional dynamics etc.)
and (ii) a precise modelling of the thermal effects themselves.
While the first problem basically needs observational inputs
(such as thermal observations of the asteroid surfaces during
fly--by space missions), the second requires theoretical efforts.
In this note we extend the precision of the thermal modelling for
a particular, but important, class of asteroid fragments.

Probably the most uncertain modelling of the thermal effects
concerns bodies which have a two--layer structure: a high--conductive
core covered with a thin layer of a low--conductive material
(``regolith''). Though there is no direct evidence of
a regolith cover for the small, meter--sized asteroid fragments
we must consider this case as one of the possibilities (e.g.
Farinella et~al. 1998). Planetary moons such as Phobos
typically do have regolith surface layers and future fly--by missions
will investigate this possibility for small near--Earth asteroids.
Theoretical considerations also support a possibility of
electrostatically sticked dust particles on cosmic bodies (e.g.
Lee 1996).

The physical reason why the previous thermal models
might not be appropriate in the case of the regolith--covered bodies
is different for the seasonal and the diurnal components of the Yarkovsky
effect.

In the seasonal effect case, the previous models (e.g. Rubincam 1998;
Farinella et~al. 1998; Vokrouhlick\'y 1999) are inappropriate due
to the assumption of homogeneity of the thermal parameters
throughout the body. The
penetration depth of the seasonal thermal wave in the low--conductive
surface regolith is of the order $l_s^1 \approx 15$~cm. This value
is to be compared with the estimated geometrical thickness of
the surface layer (possibly $\leq 1$~cm). Because the corresponding
penetration depth of the seasonal thermal wave in the high--conductive
core of the body is much larger -- e.g. $l_s^2 \approx 2$~m for a
bare--rock body -- the seasonal temperature changes affect a much
larger volume of the body than just the surface layer of
$l_s^1$ thickness (as expected from a simple thermal theory for a
homogeneous body).
Consequently, the previously obtained estimates of the seasonal
Yarkovsky effect mobility for such two--layer asteroid fragments
may be significantly incorrect and in need of revision.

Modelling the diurnal variant of the Yarkovsky effect for
bodies with very low surface conductivity presents a different problem.
Taking into account the lunar regolith thermal data (e.g.
Langseth et~al. 1973) we obtain an estimate of $38.3/\sqrt{\omega}
\, \mu{\rm m}$ for the penetration depth of the diurnal
thermal wave ($\omega$ is the rotation frequency of the body in
rad/s). Assuming a linear dependence of the rotation period with the
characteristic size of the body (e.g. Farinella et~al. 1998),
we may further express the estimated diurnal penetration
depth as a function of the fragment's size: $34.2\,\sqrt{R}\;
\mu{\rm m}$ ($R$ in meters). For meter--sized bodies of interest
in meteoritics this depth is submillimetric. Since on this scale
the surface is likely to be dusty or particulate, that is a complex
mixture of micrometric grains and voids of the same size, the
previous models
based on the homogeneity assumption may lead to wrong results.
However, the solution of this problem is much more complicated than
in the previously mentioned case of the seasonal effect.

In this note, we consider the seasonal effect only
and we remove the homogeneity assumptions of the previous
models by taking into account the two--layer structure of the
body.

\section{Theory}
Hereafter, we use the mathematical approach and, to a large extent,
also the notation introduced in the series of papers by
Vokrouhlick\'y (1998a,b; 1999). The reader is referred
details to these references for details.
\begin{figure}[t]
%\epsfxsize=8.7cm
\epsfxsize=8.2cm
\begin{center}
\mbox{\epsfbox{f1.eps}}
\end{center}
\caption[]{Interior structure of the two--layer body discussed in the
 text. \label{f1}}
\end{figure}

We consider a linearized solution of the temperature variations in a
spherical and rotating body that revolves around a radiating
source (the Sun in our case) on a circular orbit. The body of
radius $R$ is assumed to have a thin and homogeneous surface
layer of thickness $h$, characterized by thermal
conductivity $K_1$, thermal capacity $C_1$ and density $\rho_1$.
The corresponding parameters of the core are $K_2$,
$C_2$ and $\rho_2$. As noted above, we assume that the thermal
conductivities $K_1$ and $K_2$ are different; typically
$\xi_2 \equiv K_1/K_2 \simeq 1/50$. For further use we also introduce
the ratio of the thermal depths in the two parts of the body
as $\xi_1 \equiv \sqrt{K_1 \rho_2 C_2/K_2 \rho_1 C_1} \simeq
1/15$. Figure~\ref{f1} shows the geometrical structure of the
bodies discussed in this note.

Since we shall use the linearized heat conduction theory, the
temperature $T$ throughout the body is assumed to be oscillate
around a constant mean value $T_{\rm av}$: $T = T_{\rm av} + \Delta T$.
The constant average part is defined by the formal conservation of the
absorbed and reradiated energy without taking into account the
heat conduction effects. The heat conduction problem then means
to solve the temperature variation $\Delta T$ at any time $t$
everywhere in the body. The volume elements in the body are to be
parametrized by the spherical coordinates $(r;\theta,\phi)$.

A suitable scaling of the variables simplifies the mathematical
formulation of the problem. The radial coordinate $r$ is scaled
by the thermal penetration depth $l_s^1 = \sqrt{K_1/\rho_1 C_1
n}$ of the surface layer: $r \rightarrow r'=r/l_s^1$ (scaled
quantities will be denoted by prime; $n$ denotes the mean motion of the
fragment's revolution around the Sun). The temperature $\Delta T$ will
be scaled by the subsolar surface temperature $T_\star$, defined from
$\epsilon\sigma T_\star^4=\alpha {\cal E}_\star$ ($\epsilon$ is the
thermal emissivity and $\alpha$ the optical absorptivity of the
surface, $\sigma$ the Stefan--Boltzmann constant and ${\cal E}_\star$
is the solar radiation flux at the fragment's distance from the Sun).
Finally, the time $t$ will be replaced by a complex variable $\zeta=
\exp(i \lambda)$, where $\lambda = n(t-t_0)$ is the mean longitude
along the orbit. The time origin $t_0$ is
chosen so that the solar colatitude $\theta_0$ in the body's
reference system (measured from the spin axis) is such that
$\cos\theta_0=
\sin\gamma \sin\lambda$ ($\gamma$ is the fragment's rotation axis
obliquity).

The linearized heat conduction problem then satisfies the Fourier
equation (e.g. Vokrouhlick\'y 1998a; 1999)
\begin{eqnarray}
 i \zeta {\partial \over \partial \zeta}\,& &\!\!\!\!\!\!\!\!\!
  \Delta T'(r';\theta,\phi;\zeta) = \label{one} \\&=&
   {1\over r'^2}\left\{{\partial \over \partial r'}
  \left(r'^2{\partial \over \partial r'}\right) + \Lambda\left(\theta,
  \phi\right)\right\}\, \Delta T'(r';\theta,\phi;\zeta) \; , \nonumber
\end{eqnarray}
with the operator $\Lambda(\theta,\phi)$ given by
\begin{equation}
 \Lambda\left(\theta,\phi\right) = {1\over \sin\theta}\left[{\partial
  \over \partial \theta}\left(\sin\theta{\partial \over \partial \theta}
  \right)+ {1\over \sin\theta}{\partial^2 \over \partial \phi^2}\right]
 \; . \label{two}
\end{equation}
The general solution of Eq.~(\ref{one}) is constrained by three boundary
conditions. The energy flux conservation at the fragment's surface
($r'=R'$) reads
\begin{equation}
 \sqrt{2} \Delta T'(R') + \Theta\left[{\partial\Delta T'(R')\over
  \partial r'}\right] = \Delta {\cal E}'\; ,\label{three}
\end{equation}
with the incident radiation flux ${\cal E}$ normalized by the ``subsolar
value'' ${\cal E}_\star$ (${\cal E}' = {\cal E}/{\cal E}_\star$) and
$\Theta = \sqrt{K_1 \rho_1 C_1} \sqrt{n}/\epsilon\sigma T_\star^3$
defined as
the seasonal thermal parameter of the surface low--conductive layer.
A dependence of the boundary condition (\ref{three}) on the particular
surface element $(\theta,\phi)$ and time $(\zeta)$ is assumed.
The right--hand side term in Eq.~(\ref{three}) $(\Delta {\cal E}'=
{\cal E}' - \frac{1}{4})$ can be developed in the
spherical harmonics series
\begin{equation}
 \Delta {\cal E}' = \sum_{n\geq 1}\,\sum_{k=-n}^n \varepsilon_{nk}
  (\zeta)\, Y_{nk}(\theta,\phi)\; , \label{four}
\end{equation}
from which only the dipole zonal part will be relevant below.
Vokrouhlick\'y (1999) showed that the corresponding amplitude
$\varepsilon_{10}(\zeta)$ reads
\begin{equation}
 \varepsilon_{10}(\zeta) = {1\over 2i}\sqrt{{\pi \over 3}}\,
  \sin\gamma \left(\zeta-\zeta^{-1}\right)\; . \label{five}
\end{equation}
A new type of constraint is due to the boundary between
the surface low--conductive layer and the core (i.e. at $r'=R'-h'$).
Both the temperature and the heat flux are to be continuous on this
surface, thus satisfying
\begin{eqnarray}
 \Delta T'(R'-h')_+ &=&\  \Delta T'(R'-h')_- \; , \label{six} \\
 \vphantom{13mm}
 \xi_2 \left[{\partial \Delta T'(R'-h')\over
  \partial r'}\right]_+ &=&\ \left[{\partial \Delta T'(R'-h')
  \over\partial r'}\right]_- .  \label{seven}
\end{eqnarray}
The indexes $\pm$ indicate that the corresponding quantity has to be
computed
as a limit $r' \rightarrow R'-h'$ from the surface layer ($+$) and the
core ($-$) region.

Thanks to the spherical geometry the temperature inside the body may
be developed as a spherical harmonics series
\begin{equation}
 \Delta T'(r';\theta,\phi;\zeta) = \sum_{n\geq 1}\,\sum_{k=-n}^n
  t'_{nk}(r';\zeta)\,Y_{nk}(\theta,\phi)  \; , \label{eight}
\end{equation}
with the coefficients $t'_{nk}(r';\zeta)$ depending on the scaled
radial distance $r'$ from the center of the body and the time
$\zeta$. Only the dipole part $t'_{10}(r';\zeta)$ is relevant to
compute of the seasonal Yarkovsky acceleration ${\bf a}$
(e.g. Vokrouhlick\'y 1999), since
\begin{equation}
 {\bf a}(\zeta) = -{4\over 3}\sqrt{2\over 3\pi}\,\alpha\Phi\,
  t'_{10}(R';\zeta)\, {\bf s}\; . \label{nine}
\end{equation}
Here, ${\bf s}$ is the unit vector in the direction of the fragment's
spin axis and $\Phi=\pi R^2{\cal E}_\star/mc$ is the usual radiation
pressure factor ($m$ is the mass of the body and $c$ the light
velocity).

After a great deal of algebra we obtain
\begin{equation}
 t'_{10}(R';\zeta) = \sqrt{\pi\over 6}\, \sin\gamma \,
  {E_{R'} \sin(\lambda + \delta_{R'})\over 1+\chi}\; ,
  \label{ten}
\end{equation}
where
\begin{equation}
  E_{R'} \exp(i \delta_{R'})= {\beta_4 - \beta_3\over \beta_1 \beta_4
   - \beta_2 \beta_3}  \; . \label{eleven}
\end{equation}
The auxiliary complex $\beta$--factors are given by
\begin{eqnarray}
 \beta_1 &=& 1+{\chi\over 1+\chi} \psi_1(z_1)\; , \label{twelve}\\
 \beta_2 &=& 1-{3\chi\over 1+\chi}+{\chi\over 1+\chi} \psi_2(z_1)\; ,
  \label{thirteen}\\
 \beta_3 &=& \left[\xi_2-1+\xi_2\psi_1\left(z_2\right) -
  \psi_1\left(z_3\right)\right]\, \psi_3\left(z_2\right)\; ,
  \label{fourteen}\\
 \beta_4 &=& \left[-2\xi_2-1+\xi_2\psi_2\left(z_2\right) -
  \psi_1\left(z_3\right)\right]\, \psi_3\left(z_1\right)
  \; , \label{fifteen}
\end{eqnarray}
with
\begin{eqnarray}
 \psi_1(z)&=& {z\over j_1(z)} {d j_1(z)\over dz} -1\; ,
  \label{sixteen} \\
 \psi_2(z)&=& {z\over n_1(z)} {d n_1(z)\over dz} +2\; ,
  \label{seventeen} \\
 \psi_3(z)&=& {j_1(z)\over n_1(z)}\; . \label{eighteen}
\end{eqnarray}
Here, $j_1(z)$ and $n_1(z)$ are the spherical Bessel functions of the
first and the second kind. The complex arguments of the
$\beta$--functions read: $z_1 = \sqrt{-i}\,R'$, $z_2 = \sqrt{-i}\,
(R'-h')$ and $z_3 = \sqrt{-i}\, \xi_1 (R'-h')$. Finally, the real
parameter $\chi$ is defined by: $\chi=\Theta/\sqrt{2}R'$.

Using the corresponding Gauss equation we may then estimate the
secular rate of the semimajor axis evolution as:
\begin{equation}
 {da\over dt} = \frac{4\alpha}{9}\,{\Phi\over n}\,{E_{R'}\sin
  \delta_{R'}\over 1+\chi}\, \sin^2\gamma\; . \label{nineteen}
\end{equation}
This result is formally identical to Eq.~(26) in Vokrouhlick\'y
(1999) that has been derived for a homogeneous spherical fragment.
However, in our case the amplitude $E_{R'}$ and the seasonal thermal
phase lag $\delta_{R'}$ are given by the much more complex expressions
(\ref{eleven}) -- (\ref{eighteen}) above (see Appendix~A for
a numerical method that enables the evaluation of the $E_{R'}\,\sin
\delta_{R'}$ factor for $R'\gg 1$). Obviously, taking the
limit of a homogeneous body (i.e. $K_1 = K_2$, $\rho_1=\rho_2$ and
$C_1=C_2$) the semimajor axis secular drift (\ref{nineteen})
reduces to the simpler result (26) in Vokrouhlick\'y (1999).
The same applies for the limit $h\rightarrow R$, when the surface
layer formally extends over the whole body, and the limit $h
\rightarrow 0$ when the surface layer disappears.

\section{Examples and conclusions}
To assess the importance of the previously discussed effects related to
the two--layer structure of asteroid fragments, we have
performed two tests.
\begin{figure}[t]
\epsfxsize=8.7cm
\begin{center}
\mbox{\epsfbox{f2.eps}}
\end{center}
\caption[]{The average time (in Myr) necessary to move an asteroid
 fragment by $0.05$~AU in the semimajor axis vs. its radius $R$
 (in meters) by the seasonal variant of the Yarkovsky effect.
 The complete solution of this paper for three thicknesses
 of the low--conductive regolith layer (thick solid lines; $h=1$~mm,
 $h=5$~mm, $h=10$~mm) is compared with the simplified solution for a
 homogeneous body (thin solid lines; $K_1=0.0015$~W/m/K corresponding
 to the
 surface regolith -- curve~1, $K_2=1$~W/m/K corresponding to the highly
 conductive core -- curve~2). The asymptotic behavior, based on the
 solution of the 1-D heat diffusion problem (see Appendix~B)
 is shown by the dashed lines. Timescales longer than the estimated
 lifetime of the fragments are denoted by the shaded upper region.
 \label{f2}}
\end{figure}
First, we have estimated the timescale required to move a main--belt
body by $0.05$~AU due to the seasonal Yarkovsky effect, assuming a body
initially at $2.25$~AU (for a similar discussion see Hartmann
et~al. 1999). The statistical method of
Farinella \& Vokrouhlick\'y (1999) was used to include the effect
of random changes in the spin axis orientation due to collisions
with the background population of main--belt bodies. The results,
shown in Figure~\ref{f2}, indicate that the previous, simplified
approach significantly overestimates the seasonal effect mobility
of meter--sized, regolith--covered fragments, while at the same
time it significantly underestimates the seasonal effect mobility
of the dekameter--sized (and larger) bodies with the same structure.
Both results are easy to understand. We recall that the assumed
geometrical thickness of the regolith--layer (up to $10$~mm in
Figure~\ref{f2}) is smaller than the penetration depth $l_s^1$ of the
seasonal thermal wave corresponding to the thermal parameters
of the regolith (about $15$~cm). The surface temperature variations
thus penetrate well into the highly conductive core of the
body. Since the typical penetration depth $l_s^2$ of the seasonal
thermal wave in the core is of the order of $1$ -- $3$~meters,
the efficient heat conduction through the entire bulk of the
meter--sized body results in a significant decrease of its
semimajor axis mobility. In other words, the effect of the
highly conductive core helps the temperature to get
closer to an average value and thus decreases the efficiency of the
seasonal Yarkovsky effect. On the other hand, the high conductivity
of the sub--surface layer helps in increasing the ``average''
thermal parameter (indicating roughly how much the thermal
reemission lags behind the absorption of sunlight) in the
case of dekameter--sized and larger bodies. This effect then leads
to a faster seasonal Yarkovsky mobility in the semimajor axis,
as indicated in Figure~\ref{f2} (a shorter time is necessary to move
the body by $0.05$~AU).

Figure~\ref{f3} shows the average semimajor axis mobility of
asteroid fragments within their estimated collisional lifetime (compare
with Figure~3 in Farinella \& Vokrouhlick\'y 1999). The more
efficient mobility of dekameter--sized fragments covered with a
low--conductive surface layer is again confirmed [see the bump on
the thick curves in the strip (b)]. Bodies of this typical size
have both a high diurnal effect mobility and, on top of this,
a significantly increased seasonal effect mobility (as
discussed above).
\begin{figure}[t]
\epsfxsize=8.7cm
\begin{center}
\mbox{\epsfbox{f3.eps}}
\end{center}
\caption[]{Average semimajor axis drift (in AU) within the lifetime
 of an asteroid fragment vs. its radius $R$ (in kilometers). Both
 diurnal and seasonal effects included. The assumed collisional lifetime
 of $16.8\, \sqrt{R}$~Myr is taken from Farinella \& Vokrouhlick\'y 
 (1999). The complete solution of the seasonal effect (this paper)
 is shown for three thicknesses of the low--conductive regolith layer
 (thick solid lines; $h=1$~mm curve~1, $h=5$~mm curve~2, $h=10$~mm
 curve~3) and is compared with the simplified solution for a homogeneous
 body whose thermal parameters correspond to the surface layer
 (thin solid line). The diurnal variant of the Yarkovsky effect is the
 same in all cases. The region (b) of the figure, where the most
 important differences are found, correspond to
 Tunguska--like Earth impactors with typical sizes of $5$ -- $30$
 meters.
 \label{f3}}
\end{figure}

Let us now briefly comment on the possible astronomical importance of the
enhanced (``large bodies'') or inhibited (``small bodies'') seasonal
Yarkovsky effect mobility resulting from the improved model 
presented in this paper.

First, an increased
mobility in the semimajor axis of dekameter--size bodies
(Figure~\ref{f3}), provided they are endowed with a regolith
cover, could significantly enhance the flux of these bodies to the
main--belt resonances. Their preferential overabundance in the
resonances should then result in a corresponding effect among the
Earth--crossing objects. Rabinowitz's (1993, 1994) analysis indicates
a possibility for this effect, and Hartmann et~al. (1999) and
Vokrouhlick\'y \& Farinella (1999a) suggest the Yarkovsky
mobility of main--belt dekameter--sized fragments as the most
plausible explanation. In a future quantitative simulation of
these effects one should consider in detail all aspects that might
affect the semimajor axis mobility of bodies with critical sizes
($5$ -- $20$~meters), including the implications of this paper.

Second, the quantitative simulations of the cosmic ray exposure 
ages of different types of meteorites by Vokrouhlick\'y \& Farinella
(1999b) might depend on the mobility of the dekameter--sized
fragments. The current results indicate
a significant importance of secondary fragmentations for the
delivery of the parent bodies of meteorites (an idea that has been
pioneered by Wetherill 1985). Since the dekameter--sized
objects are the most likely immediate precursors of the meter--sized
meteoroids, it is very important to asses their abundance in
the main belt. In particular, one must determine whether the
dekameter--sized bodies succeed in escaping from the main belt
before being disrupted or whether they typically
suffer a catastrophic collision before reaching the resonances
(yielding thus a secondary population of meter--sized meteoroids).
Assume, for instance, a dekameter--sized, regolith--covered body
released at a $0.15$~AU distance from a resonance in the inner
region of the asteroid belt -- this is roughly the case of
asteroid Vesta. Using the results in Figure~\ref{f3} we may conclude
that the simplified approach predicts
a typical drift of about $0.1$~AU within the estimated lifetime
of a dekameter--sized body, which is not sufficient to reach the
resonances for many fragments released from Vesta. However, when
the more precise theory of this paper is used, the enhanced
mobility may protect the dekameter--sized bodies from being
disrupted on their way to the resonances, inhibiting thus to some
extend the importance of the secondary fragmentations.

A precise analysis of the orbital dynamics of the near--Earth asteroid
(1566)~Icarus represents a special problem where the results of
this paper might provide an element of the solution. According to Sitarski
(1991), Icarus' orbit determination might reveal an unexplained
secular decrease of the semimajor axis of about $7.5\times
10^{-4}$~AU/Myr. Given the fact that Icarus, which is about
$\approx 900$~m in diameter,
could have a regolith cover, the conclusions from this
paper are that one should be cautious in modelling the seasonal
Yarkovsky effect on this body (apart from the fact that Icarus'
orbit is very eccentric, which also poses a special problem
for the evaluation of the seasonal variant of the Yarkovsky
effect). Indeed, a first attempt to model the semimajor axis
mobility of Icarus' orbit (Vokrouhlick\'y, in preparation)
indicates that the classical Yarkovsky effect models fail
to predict the ``observed'' value as reported by Sitarski (1991).
Obviously, there might be other possibilities to explain this
discrepancy, that is (i) Sitarski's value is overestimated [see, e.g.,
the analysis of Yeomans (1994), who was not able to confirm
Sitarski's finding], (ii)
Icarus' size is smaller. As additional possibility we
may mention that the seasonal effect influence has been
underestimated and a fully numerical model that would simulate
consistently the abrupt change in the thermal conductivity
of the surface layer (similar to the linearized analysis in
this paper) is required. Although this possibility seems plausible
to us, we must await a detailed re--analysis of Icarus'
orbit. In any case, Icarus' dynamics might offer a unique
test of the Yarkovsky effect and thus deserves a careful
attention.

\acknowledgements{The authors thank P.~Farinella for helpful
 discussions and improvements of the final version of this paper.}
\bigskip

\noindent{\bf References}
\smallskip
\begin{description}
\item Farinella P., D. Vokrouhlick\'y, W.K. Hartmann, 1998,
 Icarus 132, 378
\item Farinella P., D. Vokrouhlick\'y, 1999, Science 283, 1507
\item Hartmann W.K., P. Farinella, D. Vokrouhlick\'y, S.J.
 Weidenschilling, A. Morbidelli, F. Marzari, D.R. Davis,
 E. Ryan, 1999, Meteoritics Planet. Sci. 34, in press
\item Langseth M.G., S.J. Keihm, J.L. Chute, 1973, Heat flow
 experiments, in: {\it Apollo 17 --- Preliminary Science Report},
 NASA SP--330
\item Lee P., 1996, Icarus 124, 181
\item Rabinowitz D.L., 1993, ApJ 407, 412
\item Rabinowitz D.L., 1994, Icarus 111, 364
\item Rubincam D.P., 1995, J. Geophys. Res. 100, 1585
\item Rubincam D.P., 1998, J. Geophys. Res. 103, 1725
\item Sitarski G., 1991, AJ 104, 1226
\item Vokrouhlick\'y D., 1998a, AA 335, 1093
\item Vokrouhlick\'y D., 1998b, AA 338, 353
\item Vokrouhlick\'y D., 1999, AA 344, 362
\item Vokrouhlick\'y D., P. Farinella, 1998, AJ 116, 2032
\item Vokrouhlick\'y D., P. Farinella, 1999a, in: {\it The Impact
 of Modern Dynamics in Astronomy}, Proceedings of the IAU Colloquium
 No.~172 (Kluwer, Dordrecht), in press
\item Vokrouhlick\'y D., P. Farinella, 1999b, to be submitted
\item Yeomans, D., 1994, AJ 104, 1266
\item Wetherill, G., 1985, Meteoritics 20, 1
\end{description}
\bigskip

\noindent{\bf Appendix A: Evaluation of the term (11)}
\medskip

\noindent
In this appendix we comment on the numerical evaluation of the
complex function $E_{R'}\,\exp(i\delta_{R'})$ in Eq.~(\ref{eleven}),
which is fundamental for estimating
the semimajor axis drift (\ref{nineteen}). Curiously, this problem
is rather tricky for reasons that will be explained below.

First, a more suitable representation of (\ref{eleven}) reads
\begin{equation}
  E_{R'} \exp(i \delta_{R'})= {1\over \beta_1+\beta_3 {\beta_1-\beta_2
   \over \beta_4- \beta_3}} \; . \label{aone}
\end{equation}
where the ratio of the $\beta$--functions in the second term of
the denominator is of particular concern. Using simple algebraic
identities for the spherical Bessel functions one proves that
\begin{equation}
 \beta_1-\beta_2 = {\chi\over 1+\chi}\, {z_1\over \psi_3(z_1)}\,
  {d\psi_3(z_1)\over dz} \; . \label{atwo}
\end{equation}
More involved algebra then leads to
\begin{eqnarray}
 \beta_4-\beta_3 &=& -z_1\, {d\psi_3(z_1)\over dz}\,
  \biggl[\xi_2 + {h\over R}\,\alpha_1\left(z_1;R,\xi_1,\xi_2\right)
  \label{athree} \\
 & & \qquad\qquad\qquad\quad\;  + \left({h\over R}\right)^2\alpha_2
  \left(z_1;R,\xi_1,\xi_2\right) + \ldots\biggr]\; , \nonumber
\end{eqnarray}
where $\alpha_1$, $\alpha_2$ etc. are some complex functions.
Notice that both the numerator $(\beta_1-\beta_2)$ and the denominator
$(\beta_4-\beta_3)$ factors are proportional to the first derivative
of the $\psi_3$ function defined in (\ref{eighteen}). One easily
proves that this term decays to zero extremely fast when the norm
$|z_1|$ of its argument increases. In quantitative terms we estimate
\begin{equation}
 {d\psi_3(z_1)\over dz} \propto \exp\left(-\sqrt{2}R'\right)
  \; . \label{afour}
\end{equation}
Since a typical value of the penetration depth of the seasonal
thermal wave in the regolith layer is $l_s^1 \approx 15$~cm and
the radius $R$ of the bodies of interest may be of the order of
$100$~meters (say), we obtain $|d\psi_3(z_1)/dz| \approx
10^{-400}$! While both the numerator and the denominator factors are
of this order, their ratio is a finite number. Obviously, without
factorizing the $d\psi_3(z_1)/dz$ term, precise evaluation
of the complex number $E_{R'}\, \exp(i\delta_{R'})$ in the
left--hand side of (\ref{aone}) would not be possible. In the
previous examples
have computed the $(h/R)$--development in (\ref{athree}) up to the
third power by using a Taylor expansion of the $\beta_3$ and $\beta_4$
functions.
\bigskip

\noindent{\bf Appendix B: Plane--parallel solution}
\medskip

\noindent As it was explained in the previous Appendix~A, evaluating
the seasonal semimajor axis drift (\ref{nineteen}) requires a
special attention for large bodies ($R' \gg 1$). In this limit,
the plane--parallel (``$1$--D'') approximation of the heat diffusion
problem might become appropriate. This approach has been developed
by Rubincam (1995) and, in the notations of this paper, by
Vokrouhlick\'y (1998b). The corresponding non--linearized theory can
be found in Vokrouhlick\'y \& Farinella (1998).

We have generalized the analysis in Vokrouhlick\'y (1998b) to account
for (i) the seasonal variant of the Yarkovsky effect, and (ii) the
composite, two--layer bodies studied in this paper. The final
formula for the semimajor axis drift due to the seasonal variant
of the Yarkovsky effect is identical with (\ref{nineteen}), provided
the $E_{R'}\,\exp(i \delta_{R'})$ function reads
\begin{equation}
 {E_{R'} \exp(i \delta_{R'})\over 1+\chi}= -{1\over 1-
  \Theta \sqrt{-{i\over 2}} {1+\xi_{12} \exp\left(2\sqrt{-i}h'
  \right) \over 1-\xi_{12} \exp\left(2\sqrt{-i}h'\right)}}\; ,
  \label{bone}
\end{equation}
keeping the same notation as above and
\begin{equation}
 \xi_{12}= {\xi_1-\xi_2\over \xi_1+\xi_2} \label{btwo}
\end{equation}
[obviously, the dependence on $R'$ disappeared from Eq.~(\ref{bone})].
Interestingly, formula (\ref{bone}) may be evaluated without
any numerical problems and can check the quality of our
approximation developed in the previous Appendix~A. Results
based on (\ref{bone}) are shown in Figure~\ref{f2} by the dashed curves.
When $R' \gg 1$, we note a fairly good match with the complete
solution (solid line), while at small sizes $R'\leq 1/\xi_1$ the
plane--parallel approximation obviously fails to predict
correct results.

\end{document}
