A&A 459, 275282 (2006)
DOI: 10.1051/00046361:20065451
Yarkovsky effect on a body with variable albedo
D. Vokrouhlický^{1,2}
1  Institute of Astronomy, Charles University,
V Holesovickách 2, 18000 Prague,
Czech Republic
2 
Department of Space Studies, Southwest Research
Institute, 1050 Walnut St., Boulder, CO 80302, USA
Received 18 April 2006 / Accepted 27 June 2006
Abstract
Context. Precise observations and long orbital arcs make determining the orbit of small nearEarth asteroids a task of ever increasing difficulty. Subtle perturbations, such as relativistic terms or radiation forces, must be taken into account.
Aims. Vokrouhlický & Milani (2000) studied the orbital effects of radiation pressure on those spherical bodies with northsouth asymmetry in taking the optical surface albedo value. However, their analysis omitted the complementary effect of recoil due to thermal radiation. Here we determine this component and study the degree of compensation between these two effects.
Methods. We use an analytic approach that requires linearization of the heatdiffusion boundary condition on the surface of the body and weakly eccentric orbits. We determine both diurnal and seasonal components of the thermal force due to albedo asymmetry, and we compute the related orbitaveraged change of the semimajor axis, eccentricity, and inclination.
Results. Within the limits of our approximation and with a plausible value of thermal parameters for small nearEarth asteroids, we show that the orbitperturbation due to anisotropic reflection in the optical waveband is largely compensated by the anisotropic reflection in the thermal waveband. The resulting effect amounts to 25% of the effect in the optical at maximum.
Key words: radiation mechanisms: thermal
Our ability to accurately track the motion of small, nearEarth asteroids
has significantly increased over the past few decades thanks to groundbased
radar and interferometric techniques and also thanks to spacecraft
visits. The observations of many asteroids extend over a longer
timespan (half a century or more for a number them) and large
telescopes allow detection of very small objects (only tens of
metres in size). This strengthens the requirements on the orbit
determination model, both by the need for a careful definition
of reference systems, and transformations between them, and also by
using a force model that is accurate enough. For instance,
relativistic terms are an unsurprising part of orbital mechanics
today (e.g., Soffel et al. 2003). Another step of extending the
orbit determination model beyond the classical frame of Newtonian
dynamics has been paved by Vokrouhlický et al. (2000) who
argued that the recoil force of thermal radiation (the Yarkovsky effect)
needs to be included in the force model for the most precise
orbits. Indeed, Chesley et al. (2003) followed that paper by
actual detection of this tiny force in the motion of (6489) Golevka,
and Vokrouhlický et al. (2005a,b) argue that several more detections
should be possible soon for both solitary asteroids and binary
systems.
Vokrouhlický & Milani (2000) extended the analysis of Vokrouhlický
et al. (2000) and considered the importance of other radiationpressure
dynamical effects in the motion of small asteroids. Among others, they
pointed out that a longterm secular change of the orbital semimajor axis
could be produced by variations in albedo on the body's surface.
Their simple model accounted for a body with a spherical shape and
northsouth difference in the albedo value. The simplest model thus
assumes the albedo reads
,
where a_{0} and a_{1} are parameters and
is colatitude
measured from the north pole on the asteroid. With that assumption,
Vokrouhlický & Milani (2000) showed that the orbital semimajor axis
a undergoes a secular change
^{}

(1) 
Here n is the mean orbital motion,
with R the radius of the body, m its mass, c the light
velocity, and
W/m^{2} is the mean solar
radiation flux at the heliocentric distance equal to the semimajor axis
of the orbit. For
a unit vector directed along the asteroids
rotation pole, we also have
with
being a unit vector in the orbital plane and perpendicular
to the pericentre direction .
With
along the
orbital angular momentum, the
form a
righthanded, orbittied frame. Vokrouhlický & Milani (2000) also
showed that this effect might raise a nonnegligible orbit
displacement in some specific cases and, perhaps, a reasonable
value of a_{1}. For instance, their example considered motion
of (1566) Icarus and a_{1}=0.01; with those parameters, the optical
reflectivity anisotropy effect has been found to displace this asteroid orbit
by as much as 15 km during its close approach in June 2015.
This value is less than the orbit uncertainty by that time, but
significantly more than the observation accuracy. Thus, on the long term,
it should become important in the orbit determination analysis.
However, a possible caveat of the analysis by Vokrouhlický & Milani (2000)
is that these authors did not consider the related modification
of the thermal component of the radiative pressure, when the albedo value
varies over the asteroid's surface (note the classical models of the
Yarkovsky effect assume a constant albedo value). In the limit of a very
small surface thermal conductivity, the thermal forces mimic those produced
by the optical radiation reflection. We might thus expect that the
effect described by Vokrouhlický & Milani (2000) would receive an
oppositebehaving counterpart in the thermal band that would just
compensate for it; note that the incoming sunlight is either reflected in
the optical band or reradiated in the thermal band according to the
local value of the surface albedo.
The resulting net orbital perturbation, combining recoil effects in optical
and thermal wavebands, could be much smaller than predicted by
Vokrouhlický & Milani (2000).
In this paper we develop a detailed linear model of the Yarkovsky forces
on a spherical body with the dipole anisotropy of the albedo coefficient given
above. By computing the secular part of the semimajor axis
perturbation, we then analyse the case of a possible compensation for the
orbital perturbation in the optical band discussed by Vokrouhlický & Milani
(2000). The less important effects on eccentricity and inclination are
analysed in the Appendix.
2 The heat diffusion problem
The analytic approach below closely follows previous works by
Vokrouhlický (1998, 1999). While using the same notation and
techniques, we basically generalise results in Vokrouhlický (1999)
by assuming the surface albedo coefficient
.
To keep it simple, however, the thermal parameters are
assumed constant. The body is assumed spherical. In the next six
sections, we give a concise stepbystep version of a
linearized heatdiffusion solution, compute the recoil force
due to the thermal radiation, and determine the major secular
orbital perturbation of the semimajor axis. On contrast to the
work of Vokrouhlický & Milani (2000), we restrict our
analysis of the orbit perturbation to the lowest eccentricity
terms. This is because the solution of the heat diffusion
problem, even in its linearized version, would be fairly
complicated for all eccentricity terms included. We thus restrict
terms in eccentricity only to the first power
^{},
which are necessary for the comparison with the optical effect
(see Eq. (1)).
Distribution of temperature T, determining the thermal state of the
body, is obtained by solution of the heat diffusion equation
(often called the Fourier equation; e.g. Landau & Lifchitz 1986)

(2) 
where
is the mean bulk density, C the specific heat capacity,
and K the thermal conductivity of the particle. The uniqueness of the
solution stems from selecting appropriate boundary conditions. In the
space domain, this means regularity of T at the centre of the particle
and energy conservation at the surface. The latter condition reads

(3) 
for each of the infinitesimal surface elements with an outward
directed unit vector .
The first term here is the
thermal energy radiated per unit of time by the particle to the
space, and the second term is the energy conducted per unit of time
inside to the particle. Thermal emissivity
is for simplicity
taken to be unity throughout this paper. The second term in (3) is
the energy conducted to the body and
on the right hand side
is the radiative flux absorbed by the surface element. The
effective boundary
condition in the time domain is set by an assumption of the periodicity
of T in one revolution of the body about the centre (Sun).
To describe the solution in detail we use spherical coordinates
with the origin r=0 at the centre of the body
and colatitude
measured from its spin axis .
The origin of the longitude
is not relevant for our work.
In order to simplify the solution we chose the following set of the
nondimensional quantities (see also Spencer et al. 1989;
Vokrouhlický 1998, 1999).
With the scaling introduced, the fundamental Eqs. (2) and (3)
now take the following form:

(7) 
and

(8) 
Here,
is the usual operator where coordinates have been replaced with
their scaled values. As a result of the scaling, the system (7)
and (8) contains a single and fundamental nondimensional
parameter

(9) 
often called the "thermal parameter''; we also note that
in its numerator is the thermal inertia of the body.
2.2 Linearized approximation
A major obstacle to an analytic solution of the heat diffusion is the
nonlinear term (
)
in the surface boundary conditions (3)
or (8). A standard procedure for handling this problem is to split
T into a suitably chosen mean value
and some small
increment :
with
becoming the solvedfor function instead of T. The difficult
quartic term in the boundary condition is then approximated using
linearization
with higher orderterms in
being neglected. The most natural
choice of
is from

(10) 
where the overbar means an average value of the corresponding quantity
over (i) the body's revolution about the centre and (ii) all surface
elements. As a result,
.
Note that only the
constant part a_{0} of the albedo value contributes to
.
With the scaled variables introduced and with the mean temperature
defined in Eq. (10), the heat diffusion problem (7)
and (8) now reads



(11) 
with the operator
given by

(12) 
and the linearized boundary condition (3) reads

(13) 
Here we denote the particle radius with R and its scaled value
and
.
If the body would have a constant albedo value (a_{1}=0),
one easily shows (e.g., Vokrouhlický 1998, 1999)
that the source
term
corresponding to an impinging radiation plane wave
from the bodyfixed direction
could be decomposed
in spherical harmonics development
with coefficients

(14) 
Here

(15) 
are auxiliary coefficients, with
,
and
are Wigner's
matrixes (e.g. Wigner 1959; in the context
of orbital dynamics
are related to the inclination
functions, see Sidlichovský 1978, 1983, 1987). In the following
work, we only need the lowestdegree
:
In the general case of a body with a dipole albedo variation and residing
on an eccentric orbit, such that the impinging solar radiation flux varies
inversely proportionally to the square of heliocentric distance d, we
still have

(21) 
but the amplitudes
get modified; we also now consider
time
to be an argument of
since the direction to the Sun
in the bodyfixed frame changes as the body rotates and revolves about
the centre. After a brief bit of algebra, one obtains
while higherdegree coefficients are not needed in this work.
In Eqs. (22)(24) we introduced

(25) 
We also recall that the orbitaveraged, mean insolation of the spherical
body is
.
This is because the orbitaveraged
value, denoted by square brackets, of
and
thus
(note we work to the first powers of
the eccentricity expansions only). To obtain an explicit dependence of
on ,
one needs to develop
using elliptic expansions. As it is
a standard tool in orbital dynamics (e.g. Brouwer & Clemence 1961), we
do not explain intermediate calculations in detail but quote only the final
results.
Given the assumed spherical geometry of the particle, it is natural
to represent
in the multipole series:

(26) 
with
denoting spherical functions. Inserting this
expansion into Eq. (11), we find that the radial and timedependent
amplitudes t'_{nk} fulfill a system of decoupled equations



(27) 
and at the surface r'=R' must satisfy

(28) 
The timedependent righthand sides
of (28)
are the sourceterm amplitudes from Eq. (21).
Assuming now a particular Fourier mode in the time development of
,
namely
(with
b integer and nonzero), Eqs. (27) and (28) admit
solution
with

(29) 
which is regular at r'=0. Here we denoted
,
,
and

(30) 
Here,
j_{n}(z) is a spherical Bessel function of the degree n of
complex argument z (e.g., Abramowitz & Stegun 1970). Note that while both
and R' depend on mean motion (or any other frequency
they would be related to),
becomes independent of it. Special
attention needs to be payed to the degreeone coefficients with n=1.
In this case, we have (see also notation in Vokrouhlický 1998, 1999)
where
.
The auxiliary functions A(x), B(x),
C(x), D(x) read
A(x) 
= 

(32) 
B(x) 
= 

(33) 
C(x) 
= 

(34) 




D(x) 
= 

(35) 




It then becomes suitable to express the complex number
in amplitudephase representation
as was used
in the last row of (31). We also keep track of the
parameter b in the argument x_{b} by denoting
.
The angle
plays the role of
a thermal lag. At the limit of zero thermal conductivity we have
and thus
and
.
With the aim of computing the thermal recoil acceleration
on the
particle, we now show that only a very limited number of amplitude terms
t'_{nk} are needed, a property which greatly simplifies the analytical solution.
In agreement with previous work, we assume thermal radiation of the particle's
surface isotropic (Lambert's law). Integrating over all contributions from
infinitesimal and oriented surface elements
,
where
,
we have (e.g. Milani et al. 1987;
Bottke et al. 2002)

(36) 
Here we denote particle's mass with m and the light velocity with c;
the minus sign indicates thermal radiation recoils on the body. Again using
linearization of the T^{4} term and introducing scaled quantities
from our solution above, we obtain

(37) 
with
the typical radiation force
factor. It is now important to recall that components of the unitary vector
can be expressed using a linear combination of spherical
function of degree 1, namely

(38) 
A complex notation again yields the more comfortable expressions
and
.
With these relations and the
orthogonality property of spherical functions, we easily find that
the thermal acceleration components
(f_{X},f_{Y},f_{Z}) in an arbitrary
frame, whose axis Z is that of the particle's spin ,
read
We conclude that only dipole coefficients t'_{1k} are needed to
compute thermal acceleration in the linearized theory (e.g.
Rubincam 1998; Vokrouhlický 1998, 1999). For that reason, we may
omit computation of other coefficients below.
The second important simplification stems from our goal of determining
the longterm effect of the thermal forces on the orbital
semimajor axis a. To the first power in eccentricity e, we have

(41) 
where
is a transverse
component of the thermal acceleration, and
is the inplane component perpendicular to the pericentre direction.
The auxiliary vector
from the right hand side of Eq. (41) reads
(again only terms up to
are retained here)

(42) 
The vector components are given here in the
inertial frame, very suitable to describing the orbital motion of the
body (recall that
is the unit vector pointing to the pericentre
of the orbit,
is a unit vector directed along the orbital
angular momentum and
). The body's
spin axis
has components
(s_{P},s_{Q},s_{k})^{T} in this frame.
We note
,
where
is the obliquity, and
we also define a complex factors
for
more compact expressions.
Secular perturbations arise as orbitaverages on the right hand sides
of Gauss perturbation equations such as (41) for the semimajor
axis. Since
only depends on the first and second powers of ,
we need to
identify the terms of the same power in
in the inertialframe
development of .
However, the solution of the heat diffusion
problem, as described above, is more natural in the bodyfixed frame.
We thus undertake the other possible approach and transform
to this later frame.
Its algebraic form becomes different for the effects due to the
Zaxis force component (often called the seasonal component)
parallel to the particle spin axis direction
and the
XYplane force components (often called the diurnal components)
in the particle's equatorial plane. We thus analyse the two cases
separately in the next two sections.
2.6.1 Seasonal component
The seasonal variant of the Yarkovsky perturbation is due to the
alongspin force component f_{Z} (Eq. (40)). Assume
in the bodyfixed frame developed in Fourier
series

(43) 
Then appropriately transforming
into the bodyfixed frame
and averaging over heliocentric motion results in the following
secular drift of the orbital semimajor axis

(44) 
where
denotes the imaginary part of a complex number z,
square brackets denote the averaging as before, and

(45) 
The amplitudes
and
could be obtained using
Eq. (31) and the appropriate Fourier terms from elliptic
expansion of
in (23). We thus straightforwardly
obtain
With those results included in Eq. (44), we finally get the
seasonalvariant contribution to the secular drift of the semimajor
axis
For convenience, we note here only terms
,
dropping
the fundamental seasonal Yarkovsky term (e.g. Vokrouhlický
1999)

(49) 
Our work suggests there is no
correction
here, which agrees with the results in Vokrouhlický &
Farinella (1999).
We also note the zero thermal conductivity limit of Eq. (48)

(50) 

Figure 1:
Residual signal strength D_{a} in percent (see Eq. (64))
as a function of spin vector components (s_{P},s_{Q}) in the orbital
plane. Left panel for obliquity
smaller than
(), right panel for obliquity larger than
(). Orbital semimajor
axis and physical parameters as of the asteroid (1566) Icarus:
AU,
W/m/K, C=800 K/kg/K,
g/cm^{3}, P=2.27 hr and
m. The best estimated
position of Icarus' spin axis orientation is shown with full circle
on the right plot (
,
;
De Angelis
1995). D=0 isoline is shown by the thick solid line; dashed
isolines are for negative values of D_{a} with an increment of
1%, such that the solid dashed isolines show values D_{a}=5%
and D_{a}=10%. 
Open with DEXTER 
2.6.2 Diurnal component
For computing the diurnal components we follow the general
approach of Vokrouhlický (1999). To impose the periodicity of the
problem, we assume that the ratio of rotation frequency to the mean motion n is an integer parameter
(it is,
however, easy to generalise the results for an arbitrary real value of
m). In this case, we highlight the Fourier series of
in the following form:

(51) 
With this notation introduced, we finally obtain the diurnal
contribution to the secular change of the semimajor axis given by



(52) 
where
denotes the real part of a complex number z, and

(53) 
After some complicated algebra Eqs. (31) and (24) yield

= 

(54) 







(55) 

= 

(56) 

= 

(57) 







(58) 

= 

(59) 
These must be included in (52) to obtain the secular
semimajor axis drift. In order to make the final result simpler,
we assume ;
i.e. the rotation frequency
is much
higher than the mean motion n (generally satisfied for asteroids
and meteoroids). In this limit we may write
so we finally have



(60) 
where again we dropped the fundamental diurnal drift rate (e.g.
Vokrouhlický 1998, 1999)

(61) 
In the limit of zero thermal conductivity, Eq. (60) becomes

(62) 
Together with the seasonal effect contribution, Eq. (50), we
thus have

(63) 
Unsurprisingly, the value of the semimajor axis drift compensates for
the perturbation (1) due to the reflection of sunlight in the
optical exactly (Vokrouhlický & Milani 2000). However, when a nonzero thermal
conductivity is taken into account, compensation of the two effects is
not exact. In the next section we investigate the magnitude of the
residual effect for conductivity values appropriate for small
asteroids.
3 Discussion and conclusions
In order to make our work quantitative, we introduce the residual signal
strength

(64) 
which gives the total albedoasymmetry effect expressed as a fraction of
the drift
due to the anisotropically
reflected sunlight in the optical as previously analysed by Vokrouhlický &
Milani (2000). As mentioned before, we have
when
,
but D_{a} is nonzero
for finite values of K; likewise, D_{a} is a function of the spin
axis direction
and
.
Following results in Vokrouhlický & Milani (2000), we chose (1566)
Icarus as an exemplary case. Figure 1 shows isolines of
D(s_{P},s_{Q}) for obliquities both smaller (left) and larger (right)
than .
We used the best estimate of the orbital and physical
parameters of this asteroid, in particular a moderately high
thermal inertia as suggested by analysis of Harris (1998). Such a
value is, however, not surprising for kilometresized nearEarth
asteroids (e.g., Delbó et al. 2006). The actual
value of s_{Q} is close to unity (symbol in the right panel of
Fig. 1), which explains the large orbital effect of
the putative albedo asymmetry predicted by Vokrouhlický & Milani
(2000). We note that D_{a} is typically negative and amounts to
16% at maximum, and for a little less than half of the phase space
of the spinaxis orientation is larger than 5%. We checked that
the maximum of D_{a} increases to nearly 25% for higher thermal
conductivity K=1 W/m/K. We also checked that D_{a} is not sensitive
to m, provided the rotation period of the asteroid is varied
within reasonable limits. For instance, a ten times longer rotation
period for Icarus would not change our result in any significant way.
If the % derived from linearized
eccentricity theory is correct, the 15 km orbit displacement of
the Icarus' orbit with respect to the Earth during the 2015 close
approach reported by Vokrouhlický & Milani (2000) might
shrink to 750 m only. If so, the suspected compensation for
the albedo anisotropy effects
in the optical and thermal bands is large enough to make the orbital
perturbation currently negligible.
Obviously, our analysis only accounts for the firstorder eccentricity
effect. In addition to large s_{Q}, the large orbital
eccentricity is also responsible for the large orbital effect of
a putative albedo variation in the Icarus case: note the (1e^{2})term in the denominator
of Eq. (1) that increases the effect by about a factor
3. While the analytic analysis in this paper
does not allow us to compute the role of highereccentricity terms
that could be determined only with a complete numerical model,
we hypothesize that the D_{a}fraction of the total orbital effect
might remain roughly the same as seen in Fig. 1.
This conjecture, however, remains to be proved by a detailed
numerical analysis with the goal of including not
only the full orbital effects but possibly also differences
in optical and thermal waveband directional reflectancies
(for the infrared band see, e.g., Lagerros 1996, 1998; or
Emery et al. 1998).
The previous example of (1566) Icarus indicates that the
orbital effect of the putative surface albedo variation is
generally negligible, about an order of magnitude smaller than
predicted by Vokrouhlický & Milani (2000). A special class
of Aten asteroids with a low value of semimajor axis and very
high eccentricity might be an exception. Due to
their proximity to the Sun and their typically small size, their surfaces
are expected to have high thermal conductivity. Interest in
this group of about 2030 asteroids was recently shown by
Margot (2003), who suggested their precise orbital tracking might
offer a valuable test of the general relativity theory. This
task strengthens requirements on the accuracy of their orbital
determination so that radiative effects should be carefully
analysed both in the optical and infrared wavebands. Their uncertainty,
for instance from limits on albedo variations over the surface
studied here, might contribute to the uncertainty budget with
which the relativistic postNewtonian parameters could be
derived from their orbit determination.
Acknowledgements
I thank Alessandro Morbidelli for helpful comments on the first
version of this paper. Part of this work was supported by
the Czech Grant Agency through grant GACR 205/05/2737.
So far, we have been dealing with the semimajor axis perturbation.
This is because, for orbits of moderately large eccentricity, this
term strongly dominates the orbital displacement, eventually observable
using precise astrometry. For orbits of lowenough eccentricity,
though, the secular perturbation of the eccentricity might produce
effects comparable to or larger than the perturbation of the semimajor
axis. In this limit, however, both effects are smaller and perhaps of
little importance. For this reason, we give just a concise overview
of the principal eccentricity and inclination effects here.
The Gauss equation for the change of osculating eccentricity reads
(e.g., Bertotti et al. 2003)

(A.1) 
where we retained the zero eccentricity terms on the right hand side and
used the notation as in the main text of the paper (additionally
defining
). Vokrouhlický & Milani (2000)
showed that the northsouth asymmetry in optical albedo results in
a secular drift of orbital eccentricity

(A.2) 
but for sake of our discussion we kept the zero eccentricity terms, while
Vokrouhlický & Milani (2000) give secular
for any eccentricity value. For reasons discussed in Sect. 1, though,
this effect should be about compensated for by the corresponding change
in the Yarkovsky orbital perturbation when
.
In what follows,
we analyse the degree of this compensation for an arbitrary surface
conductivity value by computing the zeroorder eccentricity perturbation
by the thermal force on a body with an asymmetric albedo value.
As in the case of the semimajor axis, we consider separately the
orbital effects of the spinaligned thermal force component f_{Z} and the
equatorial thermal force components (f_{X},f_{Y}) (Sect. 2.5). The relevant
perturbations are seasonal

(A.3) 
and diurnal



(A.4) 
contributions to the total
value. In these expressions we
eliminated periodic parts and only kept the secular perturbation.
The amplitudes
and
of the seasonal and
diurnal force components from Eqs. (43) and (51)
have been explicitly given in Eqs. (47), (56), and
(59). Here we need to complement them by determining
and ,
which were not necessary for the semimajor
axis effects. A straightforward computation yields

(A.5) 
and

(A.6) 
Inserting all these terms into (A.3) and (A.4), we
obtain explicit expressions



(A.7) 
for the seasonal part, and



(A.8) 
for the diurnal part in
.
We easily check that for
the total thermalforce effect reads

(A.9) 
This exactly compensates for (A.2) due to reflected sunlight
in the optical band. For ,
though, the total thermal perturbation
is not
canceled by
,
and

(A.10) 
is a measure of the residual perturbation (see Eq. (64) for the
semimajor axis). We computed
for different orientation
of the spin axis and the orbital and thermal parameters corresponding to
(1566) Icarus (Sect. 3). We found the maximum
value was about
twice as high as the maximum D_{a} value, thus reaching some 40%.
Finally, we only outline the even less important inclination
perturbation. We start from (see, e.g., Bertotti et al. 2003)

(A.11) 
given here to the first power in eccentricity, and with
(
being the argument of pericentre) and
.
Using the necessary transformations of the
orbit normal
into the bodyfixed system, we obtain



(A.12) 
for the seasonal part, and



(A.13) 
for the diurnal part in secular
.
All coefficients in these
two equations have been given above.

Abramowitz, M., & Stegun, I. A. 1970, Handbook of Mathematical
Functions (New York: Dover Publications)
(In the text)
 Bertotti,
B., Farinella, P., & Vokrouhlický, D. 2003, Physics of the
Solar System (Dordrecht: Kluwer Academic Press)
(In the text)
 Bottke, W. F.,
Vokrouhlický, D., Rubincam, D. P., & Broz, M. 2002, in
Asteroids III, ed. W. F. Bottke, A. Cellino. P. Paolicchi, & R.
P. Binzel (Tucson: University of Arizona Press), 395
(In the text)
 Brouwer, D.,
& Clemence, G. M. 1961, Methods of Celestial Mechanics (New
York: Academic Press)
(In the text)
 Chesley, S.
R., Ostro, S. J., Vokrouhlický, D., et al. 2003, Science, 302,
1739 [NASA ADS] [CrossRef] (In the text)
 De Angelis,
G. 1995, Planet. Sp. Sci., 43, 649 [NASA ADS] (In the text)
 Delbó, M.,
Dell'Oro, A., Harris, A. W., Mottola, S., & Mueller, M. 2006,
Nature, submitted
(In the text)
 Emery, J. P.,
Sprague, A. L., Witteborn, F. C., et al. 1998, Icarus, 136,
104 [NASA ADS] [CrossRef] (In the text)
 Harris, A. W.
1998, Icarus, 131, 291 [NASA ADS] [CrossRef] (In the text)
 Lagerros,
J. S. V. 1996, A&A, 310, 1011 [NASA ADS] (In the text)
 Lagerros,
J. S. V. 1998, A&A, 332, 1123 [NASA ADS] (In the text)
 Landau, L. D.,
& Lifchitz, E. M. 1986, Hydrodynamics (Moscow: Nauka)
(In the text)
 Margot, J. L.
2003, BAAS, 35(4), 1039
(In the text)
 Milani, A.,
Nobili, A. M., & Farinella, P. 1987, Nongravitational
Perturbations and Satellite Geodesy (Bristol: A. Hilger)
(In the text)
 Rubincam,
D. P. 1998, J. Geophys. Res., 103, 1725 [NASA ADS] [CrossRef] (In the text)

Sidlichovský, M. 1978, Bull. Astron. Inst. Czechosl., 29,
90 [NASA ADS] (In the text)

Sidlichovský, M. 1983, Bull. Astron. Inst. Czechosl., 34,
65 [NASA ADS] (In the text)

Sidlichovský, M. 1987, Bull. Astron. Inst. Czechosl., 38,
152 [NASA ADS] (In the text)
 Soffel, M.,
Klioner, S. A., Petit, G., et al. 2003, AJ, 126, 2687 [NASA ADS] [CrossRef] (In the text)
 Spencer, J.
R., Lebofsky, L. A., & Sykes, M. V. 1989, Icarus, 78, 337 [NASA ADS] [CrossRef] (In the text)

Vokrouhlický, D. 1998, A&A, 335, 1093 [NASA ADS] (In the text)

Vokrouhlický, D. 1999, A&A, 344, 362 [NASA ADS] (In the text)

Vokrouhlický, D., & Farinella, P. 1999, AJ, 118, 3049 [NASA ADS] [CrossRef] (In the text)

Vokrouhlický, D., & Milani, A. 2000, A&A, 362,
746 [NASA ADS] (In the text)
 Vokrouhlický, D., Milani, A.,
& Chesley, S. R. 2000, Icarus, 148, 118 [NASA ADS] [CrossRef] (In the text)

Vokrouhlický, D., Capek, D., Chesley, S. R., & Ostro, S.
J. 2005a, Icarus, 173, 166 [NASA ADS] [CrossRef] (In the text)

Vokrouhlický, D., Capek, D., Chesley, S. R., & Ostro, S.
J. 2005b, Icarus, 179, 128 [NASA ADS] [CrossRef]
 Wigner, E. P.
1959, Group Theory and its Application to the Quantum Mechanics of
Atomic Spectra (New York: Academic Press)
(In the text)
Copyright ESO 2006