A&A 459, 275-282 (2006)
DOI: 10.1051/0004-6361: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 near-Earth 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 north-south 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 heat-diffusion 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 orbit-averaged 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 near-Earth asteroids, we show that the orbit-perturbation due to anisotropic reflection in the optical waveband is largely compensated by the anisotropic reflection in the thermal waveband. The resulting effect amounts to $\sim$25% of the effect in the optical at maximum.

Key words: radiation mechanisms: thermal

1 Introduction

Our ability to accurately track the motion of small, near-Earth asteroids has significantly increased over the past few decades thanks to ground-based 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 radiation-pressure dynamical effects in the motion of small asteroids. Among others, they pointed out that a long-term 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 north-south difference in the albedo value. The simplest model thus assumes the albedo reads $A(\theta) = a_0 + a_1~\cos\theta$, where a₀ and a₁ are parameters and $\theta$ 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

$$\left[\frac{{\rm d}a}{{\rm d}t}\right]_{\rm opt} = -\frac{\Phi_a}{2n}~a_1 s_Q~ \frac{e}{1-e^2}\cdot$$ (1)

Here n is the mean orbital motion, $\Phi_a = ({\cal E}_a \pi R^2/mc)$with R the radius of the body, m its mass, c the light velocity, and ${\cal E}_a\simeq 1367/a^2$ W/m² is the mean solar radiation flux at the heliocentric distance equal to the semimajor axis of the orbit. For $\vec{s}$ a unit vector directed along the asteroids rotation pole, we also have $s_Q = \vec{s}\cdot\vec{Q}$ with $\vec{Q}$ being a unit vector in the orbital plane and perpendicular to the pericentre direction $\vec{P}$. With $\vec{k}$ along the orbital angular momentum, the $(\vec{P},\vec{Q},\vec{k})$ form a right-handed, orbit-tied frame. Vokrouhlický & Milani (2000) also showed that this effect might raise a non-negligible orbit displacement in some specific cases and, perhaps, a reasonable value of a₁. For instance, their example considered motion of (1566) Icarus and a₁=0.01; with those parameters, the optical reflectivity anisotropy effect has been found to displace this asteroid orbit by as much as $\sim$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 opposite-behaving counterpart in the thermal band that would just compensate for it; note that the incoming sunlight is either reflected in the optical band or re-radiated 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 $A(\theta) = a_0 + a_1~\cos\theta$. 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 step-by-step version of a linearized heat-diffusion 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)).

2.1 Formulation and scaling

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)

$$\rho_{\rm b} C {\partial T\over \partial t} = K~\nabla^2 T ,$$ (2)

where $\rho_{\rm b}$ 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

$$\epsilon\sigma T^4 + K \left(\vec{n}\cdot{\partial T\over \partial \vec{r}}\right) = {\cal S}$$ (3)

for each of the infinitesimal surface elements with an outward directed unit vector $\vec{n}$. 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 $\epsilon$ is for simplicity taken to be unity throughout this paper. The second term in (3) is the energy conducted to the body and ${\cal S}$ 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 $(r,\theta,\phi)$ with the origin r=0 at the centre of the body and colatitude $\theta$ measured from its spin axis $\vec{s}$. The origin of the longitude $\phi$ is not relevant for our work. In order to simplify the solution we chose the following set of the non-dimensional quantities (see also Spencer et al. 1989; Vokrouhlický 1998, 1999).

• The spatial dimension will be scaled by the thermal length
  

  $$\ell_{\rm s} = \sqrt{K\over \rho_{\rm b} C n} ,$$ (4)

  
  
  where n is the body's mean motion about the centre; in particular the radial coordinate r will have the scaled value $r' = r/\ell_{\rm s}$. Note that $\ell_{\rm s}$ is also called the penetration depth of the seasonal thermal wave since it describes the e-fold damping of the temperature in the body with frequency n.
• The time t will be replaced by a complex variable $\zeta$ reading
  

  $$\zeta = \exp(\imath n t)$$ (5)

  
  
  (here $\imath=\sqrt{-1}$ is the imaginary unit; the time origin t=0 is arbitrary and is chosen to be the pericentre passage).
• Temperature T will be scaled by an auxiliary value $T_\star$ defined by
  

  $$\epsilon \sigma T_\star^4 = {\cal S}_\star = (1-a_0)~{\cal E}_a .$$ (6)

  
  
  Recall that ${\cal E}_a$ is a solar radiation flux at the distance a of the orbital semimajor axis. Note the dipole coefficient a₁ is not included in Eq. (6). The resulting non-dimensional variable will be denoted by $T' = T/T_\star$, and similarly we define $\Delta T' = \Delta T/T_\star$.
• The energy source term ${\cal S}$ on the right hand side of (3) will be scaled by ${\cal S}_\star$, thus introducing ${\cal S}'={\cal S}/{\cal S}_\star$.

With the scaling introduced, the fundamental Eqs. (2) and (3) now take the following form:

$$\imath \zeta {\partial T'\over \partial \zeta} = \nabla'^2 T' ,$$ (7)

and

$$T'^4 + \Theta \left(\vec{n}\cdot{\partial T'\over \partial \vec{r}'}\right) = {\cal S}'.$$ (8)

Here, $\nabla'$ 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 non-dimensional parameter

$$\Theta = \frac{\sqrt{K\rho_{\rm b} C}\sqrt{n}}{\epsilon\sigma T_\star^3} ,$$ (9)

often called the "thermal parameter''; we also note that $\Gamma = \sqrt{K\rho_{\rm b} C}$ 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 non-linear term ( ${\propto} T^4$) 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 $T_{\rm av}$ and some small increment $\Delta T$: $T = T_{\rm av} + \Delta T$ with $\Delta T\ll T_{\rm av}$ becoming the solved-for function instead of T. The difficult quartic term in the boundary condition is then approximated using linearization $T^4 \approx T_{\rm av}^4 + 4T_{\rm av}^3 \Delta T + \ldots$ with higher order-terms in $\Delta T$ being neglected. The most natural choice of $T_{\rm av}$ is from

$$\epsilon \sigma T_{\rm av}^4 = \overline{\cal S} = \frac{1}{4}~ (1-a_0)~{\cal E}_a ,$$ (10)

where the over-bar 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, $T_{\rm av}= T_\star/\sqrt{2}$. Note that only the constant part a₀ of the albedo value contributes to $T_{\rm av}$.

With the scaled variables introduced and with the mean temperature defined in Eq. (10), the heat diffusion problem (7) and (8) now reads

 

$$\imath \zeta {\partial \over \partial \zeta} \Delta T'(r';\theta,... ...) + \Lambda\left(\theta, \phi\right)\right\}~ \Delta T'(r';\theta,\phi;\zeta) ,$$ (11)

with the operator $\Lambda(\theta,\phi)$ given by

$$\Lambda\left(\theta,\phi\right) = {1\over \sin\theta}\left[{\... ...1\over \sin\theta}{\partial^2 \over \partial \phi^2}\right] ,$$ (12)

and the linearized boundary condition (3) reads

$$\sqrt{2}~\Delta T'+ \Theta\left({\partial \Delta T'\over \partial r'}\right)_{R'} = \Delta {\cal S}' .$$ (13)

Here we denote the particle radius with R and its scaled value $R'=R/\ell_{\rm s}$ and $\Delta {\cal S}' = {\cal S}'- \overline{{\cal S}'}$.

2.3 Source term

If the body would have a constant albedo value (a₁=0), one easily shows (e.g., Vokrouhlický 1998, 1999) that the source term ${\cal S}'$ corresponding to an impinging radiation plane wave from the body-fixed direction $(\theta_0,\phi_0)$ could be decomposed in spherical harmonics development ${\cal S}'=\sum_{n\geq 1}~\sum_{k=-n}^n \varepsilon_{nk}(\theta_0,\phi_0)~Y_{nk}(\theta,\phi)$ with coefficients

$$\varepsilon_{nk} = {\cal J}_n~\Delta^{(n)}_{0k}\left(\theta_0\right) \exp\left(-\imath k\phi_0\right) .$$ (14)

Here

$${\cal J}_n = 2\pi~\int_0^1 {\rm d}\mu~\mu~Y_{n0}\left(\mu\right)$$ (15)

are auxiliary coefficients, with $\mu=\cos\theta$, and $\Delta^{(n)}_{0k}$ are Wigner's $(2n+1)\times (2n+1)$ matrixes (e.g. Wigner 1959; in the context of orbital dynamics $\Delta^{(n)}_{0k}$ are related to the inclination functions, see Sidlichovský 1978, 1983, 1987). In the following work, we only need the lowest-degree $\varepsilon_{nk}$:

$$\varepsilon_{00} = \frac{\sqrt{\pi}}{2} ,$$ (16)

$$\varepsilon_{10} = \sqrt{\frac{\pi}{3}}~\cos\theta_0 ,$$ (17)

$$\varepsilon_{1\pm 1} = \mp\sqrt{\frac{\pi}{6}}~\sin\theta_0~ \exp\left(\mp \imath\phi_0\right) ,$$ (18)

$$\varepsilon_{20} = \frac{\sqrt{5~\pi}}{16}~\left(3\cos^2 \theta_0 - 1\right) ,$$ (19)

$$\varepsilon_{2\pm 1} = \mp \frac{\sqrt{30~\pi}}{16}~\cos \theta_0~\sin\theta_0~\exp(\mp \imath\phi_0) .$$ (20)

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

$${\cal S}' = \sum_{n\geq 1}~\sum_{k=-n}^n \sigma_{nk}(\zeta)~Y_{nk}(\theta,\phi) ,$$ (21)

but the amplitudes $\sigma_{nk}(\zeta)$ get modified; we also now consider time $\zeta$ to be an argument of $\sigma_{nk}$ since the direction to the Sun in the body-fixed frame changes as the body rotates and revolves about the centre. After a brief bit of algebra, one obtains

   

$$\sigma_{00} \frac{\sqrt{\pi}}{2}~\left(\frac{a}{d}\right)^2 \left(1-\frac{2}{3}~\alpha_1~\cos\theta_0\right) ,$$ = (22)

$$\sigma_{10} \sqrt{\frac{\pi}{3}}~\left(\frac{a}{d}\right)^2 \left[\cos\theta_0-\frac{3}{8}~\alpha_1~\left(1+\cos^2 \theta_0\right)\right] ,$$ = (23)

$$\sigma_{1\pm 1} \mp\sqrt{\frac{\pi}{6}}~\left(\frac{a}{d} \right)^2 \left(1-\frac{3}{8}~\alpha_1~\cos\theta_0\right) ~\sin\theta_0~\exp(\mp \imath\phi_0) ,$$ = (24)

while higher-degree coefficients are not needed in this work. In Eqs. (22)-(24) we introduced

$$\alpha_1 = \frac{a_1}{1-a_0}\cdot$$ (25)

We also recall that the orbit-averaged, mean insolation of the spherical body is $\overline{{\cal S}'} = 1/4$. This is because the orbit-averaged value, denoted by square brackets, of $[\sigma_{00}]=\sqrt{\pi}/2$ and thus $[\sigma_{00}]~Y_{00}=1/4$ (note we work to the first powers of the eccentricity expansions only). To obtain an explicit dependence of $\sigma_{nk}$ on $\zeta$, one needs to develop $[a/d,\cos\theta_0, \sin\theta_0~\exp(\pm \imath \phi_0)]$ 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.

2.4 Solution of the linear heat diffusion problem

Given the assumed spherical geometry of the particle, it is natural to represent $\Delta T'$ in the multipole series:

$$\Delta T'(r';\theta,\phi;\zeta) = \sum_{n\geq 1}~\sum_{k=-n}^n t'_{nk}(r';\zeta)~Y_{nk}(\theta,\phi) ,$$ (26)

with $Y_{nk}(\theta,\phi)$ denoting spherical functions. Inserting this expansion into Eq. (11), we find that the radial- and time-dependent amplitudes t'_nk fulfill a system of decoupled equations

 

$$\imath \zeta {\partial \over \partial \zeta}~ t'_{nk}(r';\zeta) =... ...ial \over \partial r'}\right) - n\left(n+1 \right)\right\}~ t'_{nk}(r';\zeta) ,$$ (27)

and at the surface r'=R' must satisfy

$$\sqrt{2}~ t'_{nk}(R';\zeta) + \Theta\left({\partial t'_{nk}\over \partial r'}\right)_{(R';\zeta)} = \sigma_{nk}(\zeta) .$$ (28)

The time-dependent right-hand sides $\sigma_{nk}(\zeta)$ of (28) are the source-term amplitudes from Eq. (21).

Assuming now a particular Fourier mode in the time development of $\sigma_{nk}(\zeta)$, namely $\sigma_{nk}(\zeta) = \sigma_{nk}^b~\zeta^b$ (with b integer and nonzero), Eqs. (27) and (28) admit solution $t'_{nk}(r';\zeta)= t^{\prime b}_{nk}(r')~\zeta^b$ with

$$t^{\prime b}_{nk}(r') = \frac{\sigma_{nk}^b}{\sqrt{2}}~ \fra... ...}{1+\chi~ \frac{Z}{j_n(Z)}~\frac{{\rm d} j_n(Z)}{{\rm d}Z}} ,$$ (29)

which is regular at r'=0. Here we denoted $z=\sqrt{-\imath b}~r'$, $Z=\sqrt{-\imath b}~R'$, and

$$\chi=\frac{\Theta}{\sqrt{2}~R'} \cdot$$ (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 $\Theta$ and R' depend on mean motion (or any other frequency they would be related to), $\chi$ becomes independent of it. Special attention needs to be payed to the degree-one coefficients with n=1. In this case, we have (see also notation in Vokrouhlický 1998, 1999)

 

$$t^{\prime b}_{1k}(r') \frac{\sigma_{1k}^b}{\sqrt{2}}~ \frac{j_1(z)}{j_1(Z)}~\frac{1}{1+\chi~ \frac{Z}{j_1(Z)}~\frac{{\rm d} j_1(Z)}{{\rm d}Z}}$$ =

$$\frac{\sigma_{1k}^b}{\sqrt{2}~(1+\chi)}~\frac{j_1(z)} {j_1(Z)}~\frac{A(x_b)+\imath B(x_b)}{C(x_b)+\imath D(x_b)}$$ = (31)

$$\frac{\sigma_{1k}^b}{\sqrt{2}~(1+\chi)}~ \frac{j_1(z)}{j_1(Z)}~E(x_b)\exp\left[\imath \delta\left(x_b \right)\right] ,$$ =

where $x_b=\sqrt{2b}~R'$. The auxiliary functions A(x), B(x), C(x), D(x) read

    

$$-\left(x+2\right)- e^x\left[\left(x-2\right)\cos x - x \sin x\right] ,$$                           A(x) = (32)

$$-x - e^x\left[x\cos x +\left(x-2\right)\sin x\right] ,$$ B(x) = (33)

$$A(x)+ {\chi\over 1+\chi}$$ C(x) = (34)

$$\times \left\{3\left(x+2\right)+ e^x\left[ 3\left(x-2\right)\cos x+ x\left(x-3\right)\sin x\right]\right\} ,$$

$$B(x)+ {\chi\over 1+\chi}$$ D(x) = (35)

$$\times\left\{x\left(x+3\right) -e^x \left[x\left(x-3\right)\cos x- 3\left(x-2\right)\sin x\right]\right\} .$$

It then becomes suitable to express the complex number $(A+\imath B)/(C+\imath D)$in amplitude-phase representation $E\exp(\imath \delta)$ as was used in the last row of (31). We also keep track of the parameter b in the argument x_b by denoting $E_b\exp(\imath \delta_b)= E(x_b)\exp(\imath \delta(x_b))$. The angle $\delta<0$ plays the role of a thermal lag. At the limit of zero thermal conductivity we have $\chi \rightarrow 0$ and thus $E\rightarrow 1$ and $\delta\rightarrow 0$.

2.5 Thermal force computation

With the aim of computing the thermal recoil acceleration $\vec{f}$ 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 ${\rm d}\vec{S}=R^2~{\rm d}\Omega~ \vec{n}$, where ${\rm d}\Omega = \sin\theta~ {\rm d}\theta {\rm d}\phi$, we have (e.g. Milani et al. 1987; Bottke et al. 2002)

$$\vec{f}(\zeta) = -\int R^2{\rm d}\Omega~\frac{2}{3}\frac{\epsilon \sigma T^4}{mc}~\vec{n} .$$ (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⁴ term and introducing scaled quantities from our solution above, we obtain

$$\vec{f}(\zeta) = -{2\sqrt{2}\over 3\pi}~\Phi \int {\rm d}\Omega~ \Delta T'(R';\theta,\phi;\zeta)~ \vec{n} ,$$ (37)

with $\Phi = ({\cal S}_\star \pi R^2/mc)$ the typical radiation force factor. It is now important to recall that components of the unitary vector $\vec{n}$ can be expressed using a linear combination of spherical function of degree 1, namely

$$\vec{n} = \sqrt{\frac{2\pi}{3}} \left[\begin{array}{c} \phan... ...h~ (Y_{1-1}+Y_{11})\\ \sqrt{2}~ Y_{10} \end{array}\right] .$$ (38)

A complex notation again yields the more comfortable expressions $n_X\pm\imath ~n_Y = \mp~ 2\sqrt{\frac{2\pi}{3}}~Y_{1\pm 1}$ and $n_Z=2\sqrt{\frac{\pi}{3}}~Y_{10}$. 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 $\vec{s}$, read

  

$$f_X(\zeta) \pm \imath f_Y(\zeta) \mp {8\over 3\sqrt{3\pi}}~ \Phi~ t'_{1\mp 1}(R';\zeta) ,$$ = (39)

$$f_Z(\zeta) -{4\over 3}\sqrt{2\over 3\pi}~ \Phi~ t'_{10}(R';\zeta) .$$ = (40)

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.

2.6 Orbital averaging of the thermal force

The second important simplification stems from our goal of determining the long-term effect of the thermal forces on the orbital semimajor axis a. To the first power in eccentricity e, we have

$$\frac{{\rm d}a}{{\rm d}t} = \frac{2}{n}~\left(f_\tau+e~f_Q\right) ,$$ (41)

where $f_\tau = \vec{f}\cdot\hbox{\boldmath {$\tau$ }}$ is a transverse component of the thermal acceleration, and $f_Q = \vec{f}\cdot\vec{Q}$is the in-plane component perpendicular to the pericentre direction. The auxiliary vector $2~\left(\hbox{\boldmath {$\tau$ }} +e~\vec{Q} \right)$ from the right hand side of Eq. (41) reads (again only terms up to ${\propto} e^1$ are retained here)

$$2~\left(\hbox{\boldmath {$\tau$ }} +e~\vec{Q}\right) = \lef... ...antom{\imath~} \zeta^2+\zeta^{-2} \\ 0 \end{array}\right] .$$ (42)

The vector components are given here in the $(\vec{P},\vec{Q},\vec{k})$inertial frame, very suitable to describing the orbital motion of the body (recall that $\vec{P}$ is the unit vector pointing to the pericentre of the orbit, $\vec{k}$ is a unit vector directed along the orbital angular momentum and $\vec{Q}=\vec{k}\times\vec{P}$). The body's spin axis $\vec{s}$ has components (s_P,s_Q,s_k)^T in this frame. We note $s_k=\cos\gamma$, where $\gamma$ is the obliquity, and we also define a complex factors $s_\pm = s_P\pm\imath s_Q$ for more compact expressions.

Secular perturbations arise as orbit-averages on the right hand sides of Gauss perturbation equations such as (41) for the semimajor axis. Since $2~\left(\hbox{\boldmath {$\tau$ }} +e~\vec{Q} \right)$only depends on the first and second powers of $\zeta$, we need to identify the terms of the same power in $\zeta$ in the inertial-frame development of $\vec{f}$. However, the solution of the heat diffusion problem, as described above, is more natural in the body-fixed frame. We thus undertake the other possible approach and transform $2~\left(\hbox{\boldmath {$\tau$ }} +e~\vec{Q} \right)$ to this later frame. Its algebraic form becomes different for the effects due to the Z-axis force component (often called the seasonal component) parallel to the particle spin axis direction ${\bf s}$ and the XY-plane 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 along-spin force component f_Z (Eq. (40)). Assume $t'_{10}(R';\zeta)$ in the body-fixed frame developed in Fourier series

$$t'_{10}(R';\zeta) = \sum \varphi_n~\zeta^n .$$ (43)

Then appropriately transforming $2~\left(\hbox{\boldmath {$\tau$ }} +e~\vec{Q} \right)$ into the body-fixed frame and averaging over heliocentric motion results in the following secular drift of the orbital semimajor axis

$$\left[\frac{{\rm d}a}{{\rm d}t}\right]_{\rm th1} = -{8\over 3... ...t{2\over 3\pi}~ \frac{\Phi}{n}~\Im\left(s_+\varphi_+\right) ,$$ (44)

where $\Im(z)$ denotes the imaginary part of a complex number z, square brackets denote the $\zeta$-averaging as before, and

$$\varphi_+ = \varphi_1+e~\varphi_2 .$$ (45)

The amplitudes $\varphi_1$ and $\varphi_2$ could be obtained using Eq. (31) and the appropriate Fourier terms from elliptic expansion of $\sigma_{10}$ in (23). We thus straightforwardly obtain

  

$$\varphi_1 -\frac{1}{2}\sqrt{\frac{\pi}{6}}~\biggl[ s_-~\left(1-\frac{3}{4}~e\alpha_1 s_P+\frac{3}{16}~e\alpha_1 s_- \right)$$ =

$$+ \frac{3}{4}~e\alpha_1~\left(1+\sin^2\gamma\right) \biggr]~\frac{E_1\exp(\imath \delta_1)}{1+\chi} ,$$ (46)

$$\varphi_2 -\frac{3}{32}\sqrt{\frac{\pi}{6}}~\alpha_1 s_-^2~ \frac{E_2\exp(\imath \delta_2)}{1+\chi} \cdot$$ = (47)

With those results included in Eq. (44), we finally get the seasonal-variant contribution to the secular drift of the semimajor axis

 

$$\left[\frac{{\rm d}a}{{\rm d}t}\right]_{\rm th1} \frac{\Phi}{3n} \frac{e~\alpha_1}{1+\chi}~\biggl[ E_1\left(s_P\sin \delta_1+s_Q\cos\delta_1\right)$$ = (48)

$$+\frac{E_1}{4}\sin^2\gamma ~\left(s_P\sin\delta_1+ 3~ s_Q\cos\delta_1\right)$$

$$+\frac{E_2}{4}\sin^2\gamma~ \left(s_P\sin\delta_2- s_Q\cos\delta_2\right)\biggr] .$$

For convenience, we note here only terms ${\propto} \alpha_1$, dropping the fundamental seasonal Yarkovsky term (e.g. Vokrouhlický 1999)

$$\left[\frac{{\rm d}a}{{\rm d}t}\right] =\frac{4\Phi}{9n}~\frac{E_1~\sin \delta_1}{1+\chi}~\sin^2\gamma .$$ (49)

Our work suggests there is no ${\propto} e$ correction here, which agrees with the results in Vokrouhlický & Farinella (1999).

We also note the zero thermal conductivity limit of Eq. (48)

$$\left[\frac{{\rm d}a}{{\rm d}t}\right]_{{\rm th1},~K=0} =\fra... ...3n}~e~\alpha_1 s_Q~ \left(1+\frac{1}{2}~\sin^2\gamma\right) .$$ (50)

Figure 1: Residual signal strength Da in percent (see Eq. (64)) as a function of spin vector components (sP,sQ) in the orbital plane. Left panel for obliquity $\gamma$ smaller than $90^\circ$ ($s_k\geq 0$), right panel for obliquity larger than $90^\circ$ ($s_k\leq 0$). Orbital semimajor axis and physical parameters as of the asteroid (1566) Icarus: $a\sim 1.09$ AU, $K\sim 0.05$ W/m/K, C=800 K/kg/K, $\rho _{\rm b}= 2.5$ g/cm3, P=2.27 hr and $R\simeq 635$ m. The best estimated position of Icarus' spin axis orientation is shown with full circle on the right plot ( $s_P\sim -0.095$, $s_Q\sim 0.967$; De Angelis 1995). D=0 isoline is shown by the thick solid line; dashed isolines are for negative values of Da with an increment of 1%, such that the solid dashed isolines show values Da=-5% and Da=-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 $\Omega$to the mean motion n is an integer parameter $m=\Omega/n$ (it is, however, easy to generalise the results for an arbitrary real value of m). In this case, we highlight the Fourier series of $t'_{11}(R';\zeta)$in the following form:

$$t'_{11}(R';\zeta) = \zeta^m~\sum \tau_n~\zeta^n .$$ (51)

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

 

$$\left[\frac{{\rm d}a}{{\rm d}t}\right]_{\rm th2} = -{16\over 3\sq... ...ft(\cos^2\frac{\gamma}{2} ~\tau_-s_- +\sin^2\frac{\gamma}{2}~\tau_+s_+\right) ,$$ (52)

where $\Re(z)$ denotes the real part of a complex number z, and

$$\tau_\pm = \tau_{\pm 1}+e~\tau_{\pm 2} .$$ (53)

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

    

$$\tau_1 \frac{\imath}{2}\sqrt{\frac{\pi}{3}}~\frac{E_{m+1} \exp(\imath\delta_{m+1})}{\sin\gamma~(1+\chi)}$$ = (54)

$$\times \biggl\{s_-~\biggl[\sin^2\frac{\gamma}{2}+\frac{3}{16} ~e\alpha_1~\left(s_P\cos\gamma+\imath s_Q\right)\biggr]$$

$$-\frac{3}{16}~e\alpha_1~\sin^2\gamma~\left(\cos\gamma+ \cos^2\frac{\gamma}{2}\right)\biggr\} ,$$ (55)

$$\tau_2 \frac{3~\imath}{32}\sqrt{\frac{\pi}{3}}~\frac{E_{m+2} \exp(\imath\delta_{m+2})}{\sin\gamma~(1+\chi)}~\alpha_1 \sin^2\frac{\gamma}{2}~s_-^2 ,$$ = (56)

$$\tau_{-1} \frac{\imath}{2}\sqrt{\frac{\pi}{3}}~\frac{E_{m-1} \exp(\imath\delta_{m-1})}{\sin\gamma~(1+\chi)}$$ = (57)

$$\times \biggl\{s_+~\biggl[\cos^2\frac{\gamma}{2}-\frac{3}{16} ~e\alpha_1~\left(s_P\cos\gamma+\imath s_Q\right)\biggr]$$

$$+\frac{3}{16}~e\alpha_1~\sin^2\gamma~\left(\cos\gamma- \cos^2\frac{\gamma}{2}\right)\biggr\} ,$$ (58)

$$\tau_{-2} -\frac{3~\imath}{32}\sqrt{\frac{\pi}{3}}~\frac{E_{m-2} \exp(\imath\delta_{m-2})}{\sin\gamma~(1+\chi)}~\alpha_1 \cos^2\frac{\gamma}{2}~s_+^2 .$$ = (59)

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

 

$$\left[\frac{{\rm d}a}{{\rm d}t}\right]_{\rm th2} = -\frac{\Phi}{6... ...os\gamma}{1+\chi}~E_m~\left(s_P\sin\delta_m-s_Q\cos \gamma\cos\delta_m\right) ,$$ (60)

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

$$\left[\frac{{\rm d}a}{{\rm d}t}\right] =-\frac{8\Phi}{9n}~\frac{E_m~\sin \delta_m}{1+\chi}~ \cos\gamma .$$ (61)

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

$$\left[\frac{{\rm d}a}{{\rm d}t}\right]_{{\rm th2},~K=0} = \frac{\Phi}{6n}~e~ \alpha_1 s_Q~\cos^2\gamma .$$ (62)

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

$$\left[\frac{{\rm d}a}{{\rm d}t}\right]_{{\rm th},~K=0} = \lef... ...t}\right]_{{\rm th2},~K=0} = \frac{\Phi}{2n}~e~\alpha_1s_Q .$$ (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 non-zero 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

$$D_a = \frac{\left[{\rm d}a/{\rm d}t\right]_{\rm th1} + \left[... ...\right]_{\rm opt}}{\left[{\rm d}a/{\rm d}t\right]_{\rm opt}} ,$$ (64)

which gives the total albedo-asymmetry effect expressed as a fraction of the drift $\left[{\rm d}a/{\rm d}t\right]_{\rm opt}$ due to the anisotropically reflected sunlight in the optical as previously analysed by Vokrouhlický & Milani (2000). As mentioned before, we have $D_a\rightarrow 0$when $K\rightarrow 0$, but D_a is non-zero for finite values of K; likewise, D_a is a function of the spin axis direction ${\bf s}$ and $m=\Omega/n$.

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 $90^\circ$. 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 kilometre-sized near-Earth 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 spin-axis 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 $D\sim 5$% derived from linearized eccentricity theory is correct, the $\sim$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 $\sim$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 first-order 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 (1-e²)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 higher-eccentricity 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 20-30 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 post-Newtonian 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.

Appendix A: Secular drifts of eccentricity and inclination

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 low-enough 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)

$$\frac{{\rm d}e}{{\rm d}t} = \frac{1}{2~na}\left[3f_Q -\Im\left(\zeta^2 f_-\right) \right] ,$$ (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 $f_\pm = f_P \pm\imath f_Q$). Vokrouhlický & Milani (2000) showed that the north-south asymmetry in optical albedo results in a secular drift of orbital eccentricity

$$\left[\frac{{\rm d}e}{{\rm d}t}\right]_{\rm opt} = -\frac{\Phi_a}{3~na}~a_1 s_Q~ ,$$ (A.2)

but for sake of our discussion we kept the zero eccentricity terms, while Vokrouhlický & Milani (2000) give secular $\left[{\rm d}e/{\rm d}t\right]_{\rm opt}$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 $K\rightarrow 0$. In what follows, we analyse the degree of this compensation for an arbitrary surface conductivity value by computing the zero-order 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 spin-aligned thermal force component f_Z and the equatorial thermal force components (f_X,f_Y) (Sect. 2.5). The relevant perturbations are seasonal

$$\left[\frac{{\rm d}e}{{\rm d}t}\right]_{\rm th1} = -\frac{2}... ...} \left[3~s_Q~\varphi_0 + \Im\left(s_+\varphi_2\right)\right]$$ (A.3)

and diurnal

 

$$\left[\frac{{\rm d}e}{{\rm d}t}\right]_{\rm th2} = -\frac{4}{3\sq... ...os^2\frac{\gamma}{2}~\tau_{-2}s_- +\sin^2\frac{\gamma}{2}~ \tau_{+2} s_+\biggr]$$ (A.4)

contributions to the total ${\rm d}e/{\rm d}t$ value. In these expressions we eliminated periodic parts and only kept the secular perturbation. The amplitudes $\varphi_2$ and $\tau_{\pm 2}$ 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 $\varphi_0$ and $\tau_0$, which were not necessary for the semimajor axis effects. A straightforward computation yields

$$\varphi_0 = -\frac{3}{16}\sqrt{\frac{\pi}{6}}~\alpha_1\frac{2+\sin^2 \gamma}{1+\chi}$$ (A.5)

and

$$\tau_0 = -\frac{3}{64}\sqrt{\frac{\pi}{3}}~\alpha_1\sin 2\gamma~ \frac{\imath E_m \exp(\imath \delta_m)}{1+\chi} \cdot$$ (A.6)

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

 

$$\left[\frac{{\rm d}e}{{\rm d}t}\right]_{\rm th1} = \frac{\Phi}{8~... ...frac{1}{6} E_2\sin^2\gamma\left(s_Q \cos\delta_2-s_P\sin\delta_2\right)\Bigr] ,$$ (A.7)

for the seasonal part, and

 

$$\left[\frac{{\rm d}e}{{\rm d}t}\right]_{\rm th2} = -\frac{\Phi}{1... ...ta_m\cos\gamma -s_Q\cos\delta_m\left(1-\frac{5}{4}\sin^2 \gamma\right)\biggr] ,$$ (A.8)

for the diurnal part in ${\rm d}e/{\rm d}t$. We easily check that for $K\rightarrow 0$ the total thermal-force effect reads

$$\left[\frac{{\rm d}e}{{\rm d}t}\right]_{{\rm th}, K=0} = \frac{\Phi}{3~na}~\alpha_1 s_Q .$$ (A.9)

This exactly compensates for (A.2) due to reflected sunlight in the optical band. For $K\neq 0$, though, the total thermal perturbation $\left[{\rm d}e/{\rm d}t\right]_{\rm th1} + \left[{\rm d}e/{\rm d}t\right]_{\rm th2}$ is not canceled by $\left[{\rm d}e/{\rm d}t\right]_{\rm opt}$, and

$$D_{\rm e} = \frac{\left[{\rm d}e/{\rm d}t\right]_{\rm th1} + ... ...}t\right]_{\rm opt}}{\left[{\rm d}e/{\rm d}t\right]_{\rm opt}}$$ (A.10)

is a measure of the residual perturbation (see Eq. (64) for the semimajor axis). We computed $D_{\rm e}$ for different orientation of the spin axis and the orbital and thermal parameters corresponding to (1566) Icarus (Sect. 3). We found the maximum $D_{\rm e}$ 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)

$$\frac{{\rm d}I}{{\rm d}t} = \frac{f_k}{na}~\Re\left[\xi~\left(\zeta+\frac{e}{2} \zeta^2-\frac{3e}{2}\right)\right] ,$$ (A.11)

given here to the first power in eccentricity, and with $\xi=\exp(\imath \omega)$ ($\omega$ being the argument of pericentre) and $f_k = \vec{f}\cdot\vec{k}$. Using the necessary transformations of the orbit normal $\vec{k}$ into the body-fixed system, we obtain

 

$$\left[\frac{{\rm d}I}{{\rm d}t}\right]_{\rm th1} = -\frac{4}{3}\s... ...{-1}\left(\varphi_1+\frac{e}{2}~\varphi_2-\frac{3e}{2} ~\varphi_0\right)\right]$$ (A.12)

for the seasonal part, and

 

$$\left[\frac{{\rm d}I}{{\rm d}t}\right]_{\rm th2} = \frac{4}{3\sqr... ...t) + \xi^{-1}\left(\tau_1+ \frac{e}{2}~\tau_2-\frac{3e}{2}~\tau_0\right)\biggr]$$ (A.13)

for the diurnal part in secular ${\rm d}I/{\rm d}t$. All coefficients in these two equations have been given above.

References

• 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, Non-gravitational 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)