A&A 437, 753-764 (2005)
DOI: 10.1051/0004-6361:20053031

Synchronous motion in the Kinoshita problem

Application to satellites and binary asteroids

S. Breiter¹ - B. Melendo² - P. Bartczak¹ - I. Wytrzyszczak¹

1 - Astronomical Observatory of A. Mickiewicz University, Soneczna 36, 60-286 Poznan, Poland
2 - Grupo de Mecánica Espacial, Universidad de Zaragoza, 50009 Zaragoza, Spain

Received 9 March 2005 / Accepted 18 March 2005

Abstract
A Lie-Poisson integrator with Wisdom-Holman type splitting is constructed for the problem of a rigid body and a sphere (the Kinoshita problem). The algorithm propagates not only the position, momentum and angular momentum vector of the system, but also the tangent vector of "infinitesimal displacements''. The latter allow to evaluate the maximum Lyapunov exponent or the MEGNO indicator of Cincotta and Simó. Three exemplary cases are studied: the motion of Hyperion, a fictitious binary asteroid with Hyperion as one of the components, and the binary asteroid 90 Antiope. In all cases the attitude instability of the rotation state with spin vector normal to an equatorial orbit influences stability of the system at lower rotation rates. The MEGNO maps with variations restricted to the orbital plane for position and momentum, and to the orbit normal direction for the angular momentum resemble usual Poincaré sections. But if no restriction is imposed on the variations, some stable zones turn into highly chaotic regions, often retaining the shape of their boundaries.

Key words: chaos - methods: numerical - celestial mechanics - planets and satellites: individual: Hyperion - minor planets, asteroids

1 Introduction

The dynamics of a system consisting of a triaxial rigid body and a homogeneous sphere has a particular importance for celestial mechanics. The problem is nonintegrable in general, but still much simpler than the general problem of two rigid bodies and may serve as the first approximation when one of the bodies is almost spherical. Probably the most comprehensive source about the general problem of a rigid body and a sphere in the absence of resonances is still the work due to Kinoshita (1972). The author was obviously not the first one to discuss the problem, but he provided a first order analytical theory of orbital and rotational motion in the system. He also demonstrated that - without any approximation - an appropriate choice of the reference frame leads to the constant separation ($180^\circ$) between the ascending nodes of orbit and equator. We hope that naming the "triaxial body + sphere'' problem after H. Kinoshita is a justified abbreviation to be used in this paper and it should not be taken as a disrespect to other authors who worked on this problem.

There are two satellite-type limit cases of the Kinoshita problem: when the mass of the sphere ($m_{\rm s}$) is negligibly small, the system represents the motion of a point mass satellite around a freely rotating planet with the mass $m_{\rm e} \gg m_{\rm s}$, but in the opposite case, when $m_{\rm e} \ll m_{\rm s}$, one obtains the problem of a rotating satellite of a spherical planet. Both problems are quite common in the natural and artificial satellites studies and were extensively studied, but the discovery of binary asteroids calls for more interest in studying the intermediate cases that fall half way between the two extremes. Indeed, a significant part of binary systems, like 90 Antiope or 617 Patroclus, have components with comparable masses (Merline et al. 2002, Table 2).

Even in the triaxial satellite limit case, however, one should observe that most of the studies were performed for the planar case, with the triaxial body rotating around the axis of the maximum moment of inertia normal to the orbital plane. This case is perfectly suited for the application of Poincaré surface of section method. On the other hand, Wisdom et al. (1984) pointed out the attitude instability phenomenon that makes such restriction unjustified in some situations. One of the objectives of the present paper is to compare the results obtained with and without the planar restriction.

The present paper actually describes our research on a special case of the Kinoshita problem: the triaxial body is an ellipsoid. The ellipsoid as a source of the gravity field is approximated by a model of two material segments and a central mass - a modification of the two segments model introduced by Bartczak & Breiter (2003). Sections 2 and 3 present the formalism and a Lie-Poisson numerical integration algorithm for the Kinoshita problem. The algorithm has many common points with the work of Touma & Wisdom (1994) hence we do not repeat their discussion of the conservation of Lie-Poisson structure and of approximate conservation of the Hamiltonian. We only recall that Lie-Poisson integrators integrate Hamiltonian equations of the Lie-Poisson type in a same manner that symplectic integrators numerically solve canonical Hamiltonian equations. In Sect. 4 we extend the integrator, providing the algorithm to propagate the associated tangent vector. For the Keplerian motion we incorporate the formulation of Mikkola & Innanen (1999) and we add the recipe for the rotational component of the problem. Tangent maps allow the application of the MEGNO (Mean Exponential Growth factor of Nearby Orbits) stability indicator invented by Cincotta & Simó (2000). Section 5 recalls basic facts about MEGNO and its relation to the maximum Lyapunov exponent. MEGNO is a reliable tool originally developed for the galactic dynamics studies (Cincotta et al. 2003; Cincotta & Simó 2000) but then successfully transplanted to the dynamics of planetary systems (Gozdziewski et al. 2001) and artificial satellite problem (Breiter et al. 2005). Gozdziewski (2003) originated combining the MEGNO computation with symplectic integration of equations of motion. In the rotational dynamics problems, the MEGNO indicator was successfully applied by Pavlov & Maciejewski (2003) who studied the 3:2 spin-orbit coupling of Mercury.

Finally, Sect. 6 presents the MEGNO analysis for three sample problems. First we revisit the classical problem of the Saturn's moon Hyperion. This is an example of the Kinoshita problem in the triaxial satellite limit. Two following examples treat the comparable masses case of the Kinoshita problem intended to shed light on the dynamics of binary asteroids. The rigid body and a sphere model was already used in this context as a first approximation (Scheeres 2002). In order to provide a direct comparison with the satellite limit, we first study a fictitious system composed of Hyperion and a spherical companion whose mass is only 5 times greater than Hyperion. Then we focus on a real problem of the binary system 90 Antiope.

   
2 Equations of motion

We study the evolution of 9 variables

$$\vec{\zeta} = \left( R_1, R_2, R_3,P_1, P_2, P_3, G_1, G_2, G_3 \right)^{\rm T},$$ (1)

describing the motion of a spherical mass $m_{\rm s}$ with respect to a homogeneous ellipsoid with the mass $m_{\rm e}$. The ellipsoid has the matrix of inertia

$${\bf I} = \left(% \begin{array}{ccc} A & 0 & 0 \\ 0 & B... ...0 & a^2+c^2 & 0 \\ 0 & 0 & a^2+b^2 \\ \end{array}% \right),$$ (2)

and its momentum vector $\vec{G}$ is expressed in the body frame. The origin O of the body frame is fixed at the ellipsoid's centre of mass, whereas the orthogonal axes OX₁, OX₂ and OX₃ are directed towards the smallest, the medium and the maximum inertia axes, forming a right-handed system. Position $\vec{R}$and orbital momentum $\vec{P}$ of the sphere are also reckoned in the body frame.

The reduced mass m is a function of the ellipsoid mass $m_{\rm e}$ and the sphere mass $m_{\rm s}$

$$m = \frac{m_{\rm e}~m_{\rm s}}{m_{\rm e}+m_{\rm s}}\cdot$$ (3)

Using an appropriate Lie-Poisson bracket (Maciejewski 1996; Olver 1993), equations of motion can be derived from the Hamiltonian function

$$\mathcal{H}(\vec{\zeta}) = \frac{\vec{P}^2}{2~m} + \frac{1}{2}~\vec{G}^{\rm T} {\bf I}^{-1} \vec{G} + V(\vec{R}),$$ (4)

and read

$$\dot{\vec{\zeta}} = {\bf J} \frac{\partial \mathcal{H}}{\partial \vec{\zeta}}\cdot$$ (5)

The Lie-Poisson structure matrix ${\bf J}$ is

$${\bf J} = \left(% \begin{array}{ccc} {\bf0} & {\bf E} & {\... ... {\bf Q}(\vec{P}) & {\bf Q}(\vec{G}) \\ \end{array}% \right),$$ (6)

where ${\bf0}$ and ${\bf E}$ are $3 \times 3$ zero and unit matrices respectively. The vector product matrix ${\bf Q}$ is defined as

$${\bf Q}(\vec{F}) = \left(% \begin{array}{ccc} 0 & -F_3 & F... ... F_3 & 0 & -F_1 \\ -F_2 & F_1 & 0 \\ \end{array}% \right),$$ (7)

for any vector $\vec{F}$. Using ${\bf Q}$ we can express the vector product as a matrix product

$$\vec{A} \times \vec{B} = {\bf Q}(\vec{A})~\vec{B}.$$

With the Lie-Poisson bracket $\{f,g\}$of two arbitrary functions f,g

$$\{f,g\} = \left(\frac{\partial f}{\partial \vec{\zeta}}\right)^{\rm T} {\bf J} ~ \frac{\partial g}{\partial \vec{\zeta}},$$ (8)

equations of motion (5) become simply

$$\dot{ \vec{\zeta} } = \{ \vec{\zeta},\mathcal{H}\}.$$ (9)

A similar formalism has been used by Touma & Wisdom (1994). Our first objective is to derive a numerical integration algorithm that conserves the Lie-Poisson bracket (8).

   
3 Integration algorithm

In order to integrate numerically Eq. (5), we construct the algorithm based on a splitting method, using the Wisdom-Holman strategy of partitioning Hamiltonian in a way that introduces a small parameter as a multiplier of the truncation error (Wisdom & Holman 1991). There are two small parameters in our problem: $\varepsilon_1 = O(C-A)$ and $\varepsilon_2 = O(C-B)$, both related to the deviation of the ellipsoid shape from the spherical symmetry. Their magnitudes depend on the problem at hand, but $b \leqslant a$ implies $\varepsilon_2 \leqslant \varepsilon_1$. Thus, the potential V of the ellipsoid can be considered as the sum of a purely Keplerian part

$$V_0(R) = - \frac{k^2~m_{\rm e}~m_{\rm s}}{R },$$ (10)

and of a perturbation $V_1 = V - V_0(R) = O(\varepsilon_1)$.

The Hamiltonian (4) can be split into four separately integrable parts

$$\mathcal{H}= \mathcal{H}_0 + \mathcal{H}_1(\varepsilon_1) + \mathcal{H}_2(\varepsilon_2) + \mathcal{H}_3(\varepsilon_1),$$ (11)

where

$$\mathcal{H}_0 \frac{\vec{P}^2}{2~m} - \frac{k^2~m_{\rm e}~m_{\rm s}}{R } + \frac{\vec{G}^{\rm T} \vec{G}}{2~C},$$ = (12)

$$\mathcal{H}_1 \left( \frac{1}{A}-\frac{1}{C} \right)~\frac{G_1^2}{2},$$ = (13)

$$\mathcal{H}_2 \left( \frac{1}{B}-\frac{1}{C} \right)~\frac{G_2^2}{2},$$ = (14)

$$\mathcal{H}_3 V(\vec{R}) + \frac{k^2~m_{\rm e}~m_{\rm s}}{R } ~=~ V_1(\vec{R}).$$ = (15)

Let us introduce the Lie derivatives

$$L \vec{\zeta} = {\bf J} \frac{\partial \mathcal{H}}{\partial ... ...\partial \mathcal{H}_j}{\partial \vec{\zeta}}, ~~j=0,\ldots,3,$$

and associated exponential operators $\exp[h~L]$, $\exp[h~L_j]$, that stand for the mapping from the initial conditions $\vec{\zeta}_0 = \vec{\zeta}(t_0)$ to $\vec{\zeta}_1 = \vec{\zeta}(t_0+h)$. Accordingly, $\exp[h~L]$ describes the motion in the full problem (5)

$$\vec{\zeta}_1 = \exp[h~L]~\vec{\zeta}_0, \qquad$$ (16)

and $\exp[h~L_j]$ defines the motion under the action of the $\mathcal{H}_j$ part of the Hamiltonian.

The main building block of the Yoshida-type integrators (Yoshida 1993) is a second order "leapfrog'' map $\Psi_h$. If $\mathcal{H}= \mathcal{H}_0 + \mathcal{H}_1(\varepsilon)$, it is defined as

$$\Psi_h \vec{\zeta}_0 \exp\left[\frac{h}{2}~L_0\right] \circ \exp\left[h~L_1\right] \circ \exp\left[\frac{h}{2}~L_0\right] ~ \vec{\zeta}_0$$ =

$$\simeq \exp\left[h~L\right]~\vec{\zeta}_0 + O\left(\varepsilon~h^3\right).$$ (17)

In our case, with the Hamiltonian (11), the leapfrog map has a different form $\Phi_h$; it becomes

 

$$\Phi_h \vec{\zeta}_0 \exp\left[\frac{h}{2}~L_0\right] \circ \exp\left[\frac{h}{2}~L_1\right] \circ \exp\left[\frac{h}{2}~L_2\right] \circ \exp\left[h~L_3\right]$$ =

$$\circ \exp\left[\frac{h}{2}~L_2\right] \circ \exp\left[\frac{h}{2}~L_1\right] \circ \exp\left[\frac{h}{2}~L_0\right] ~ \vec{\zeta}_0$$

$$\simeq \exp[h~L]~\vec{\zeta}_0 + O\left(\varepsilon_1 ~h^3\right).$$ (18)

The fourth order integrator (SI4) is obtained as a composition

$$\vec{\zeta}_1 \approx \Phi_{c_1h} \circ \Phi_{c_0h} \circ \Phi_{c_1h} ~\vec{\zeta}_0 + O\left(\varepsilon_1~h^5\right),$$ (19)

where $c_0 = 1/(1- \sqrt[3]{4}~)$ and c₁= (1-c₀)/2.

In order to create the Lie-Poisson integrator, we have to find the explicit form of the maps $\exp[(h/2)~L_j]$.

3.1 Motion due to $\mathcal{H}{_0}$

  $\mathcal{H}{_0}$ describes the relative two-body problem in the reference frame rotating with a constant angular rate around the axis determined by $\vec{G}$. One may easily check that $\{\vec{G},\mathcal{H}_0\}=\vec{0}$and $\vec{G}$ is a constant vector. Accordingly, we will only discuss the motion in $\vec{R}$ and $\vec{P}$.

As a matter of fact, we can further split $\mathcal{H}{_0}$ into

$$\mathcal{H}_0 = \mathcal{A} + \mathcal{B},$$

$$\mathcal{A} = \frac{\vec{P}^2}{2~m} - \frac{k^2~m_{\rm e}~m_{\rm s}}{R},$$ (20)

$$\mathcal{B} = \frac{\vec{G}^2}{2~C},$$ (21)

and then we observe that $\{\mathcal{A},\mathcal{B}\}=0$. The two parts commute and so the motion due to $\mathcal{H}{_0}$ can be represented as pure Keplerian motion in $\vec{R}, \vec{P}$ generated by  $\mathcal{A}$

$$\dot{\vec{R}} \{\vec{R},\mathcal{A}\} ~=~ \frac{\vec{P}}{m},$$ = (22)

$$\dot{\vec{P}} \{\vec{P},\mathcal{A}\} ~=~ - \frac{k^2~m_{\rm e}~m_{\rm s}}{R^3}~\vec{R},$$ = (23)

followed by the $\mathcal{B}$ generated uniform rotation around a fixed axis

$$\dot{\vec{R}} \{\vec{R},\mathcal{B}\} ~=~ \frac{\vec{R} \times \vec{G}}{C},$$ = (24)

$$\dot{\vec{P}} \{\vec{P},\mathcal{B}\} ~=~ \frac{\vec{P} \times \vec{G}}{C}\cdot$$ = (25)

As it follows, $\exp[\tau L_0] = \Phi^{{\rm rot}}_\tau \circ \Phi^{{\rm Kep}}_\tau$ and the sequence of the two maps is not important.

The Keplerian motion $\Phi^{{\rm Kep}}_\tau$ can be most easily implemented in terms of the Gaussian functions $f(\tau)$ and $g(\tau)$(Mikkola & Innanen 1999)

$$\vec{R}_1 f~\vec{R}_0 + \frac{g}{m} ~ \vec{P}_0,$$ = (26)

$$\vec{P}_1 m~\dot{f}~\vec{R}_0 + \dot{g}~\vec{P}_0.$$ = (27)

The rotation $\Phi^{{\rm rot}}_\tau$ should be performed around the axis specified by the unit vector $\vec{\hat{G}} = G^{-1} \vec{G}$, with the rotation angle $\theta = G~C^{-1}~\tau$. This can be achieved by means of the Rodrigues formula (Marsden & Ratiu 1999)

$$\vec{R}_1 = {\bf M}_0 \vec{R}_0, \qquad \vec{P}_1 = {\bf M}_0 \vec{P}_0,$$ (28)

where

$${\bf M}_0 = {\bf E} + \sin{\theta} ~ {\bf Q}(\vec{\hat{G}}) + 2 \sin^2{(\theta/2)} ~ {\bf Q}^2(\vec{\hat{G}}).$$ (29)

3.2 Motion due to $\mathcal{H}{_1}$ and $\mathcal{H}{_2}$

The Hamiltonian function $\mathcal{H}{_1}$ generates the rotation around the axis OX₁ of the body frame

$$\dot{\vec{\zeta}} = \{ \vec{\zeta}, \mathcal{H}_1\} = G_1~... ...{X}_1 \\ \vec{G} \times \vec{X}_1 \\ \end{array}% \right) ,$$ (30)

where $\vec{X}_1 = (1,0,0)^{\rm T}$. Thus, the map $\exp[\tau L_1]$ amounts to

$$\vec{R}_1 = {\bf M}_1(\varphi_1) \vec{R}_0, ~~~ \vec{P}_1 = {... ..._1) \vec{P}_0, ~~~ \vec{G}_1 = {\bf M}_1(\varphi_1) \vec{G}_0,$$ (31)

where ${\bf M}_1(\varphi_1)$ is the elementary rotation matrix

$${\bf M}_1(\varphi_1) = \left(% \begin{array}{ccc} 1 & 0 & ... ...& -\sin{\varphi_1} & \cos{\varphi_1} \\ \end{array}% \right),$$ (32)

and the rotation angle is

$$\varphi_1 = G_1~\frac{C-A}{A~C}~\tau.$$ (33)

Similarly to the previous case, $\mathcal{H}{_2}$ generates the rotation around the axis OX₂

$$\dot{\vec{\zeta}} = \{ \vec{\zeta}, \mathcal{H}_2\} = G_2~... ...{X}_2 \\ \vec{G} \times \vec{X}_2 \\ \end{array}% \right) ,$$ (34)

and $\exp[\tau L_2]$ is defined as

$$\vec{R}_1 = {\bf M}_2(\varphi_2) \vec{R}_0, ~~~ \vec{P}_1 = {... ..._2) \vec{P}_0, ~~~ \vec{G}_1 = {\bf M}_2(\varphi_2) \vec{G}_0,$$ (35)

where

$${\bf M}_2(\varphi_2) = \left(% \begin{array}{ccc} \cos{\va... ...sin{\varphi_2} & 0 & \cos{\varphi_2} \\ \end{array}% \right),$$ (36)

with

$$\varphi_2 = G_2~\frac{C-B}{B~C}~\tau.$$ (37)

3.3 Motion due to $\mathcal{H}{_3}$

The perturbing part of the ellipsoid's potential $\mathcal{H}{_3}$depends on the $\vec{R}$ vector only and it leads to trivial equations

$$\dot{\vec{\zeta}} = \{ \vec{\zeta}, \mathcal{H}_3\} = \left... ...{R} \times \nabla V_1 \\ \end{array}% \right) = {\rm const.},$$ (38)

where $\nabla V_1 = \frac{\partial V_1}{\partial \vec{R}}$.

In these circumstances we obtain the simple "kick map'' $\exp[\tau L_3]$

 

$$\vec{R}_1 \vec{R}_0$$ =

$$\vec{P}_1 \vec{P}_0 - \nabla V_1(\vec{R}_0)~\tau,$$ = (39)

$$\vec{G}_1 \vec{G}_0 + \left(\vec{R}_0 \times \nabla V_1(\vec{R}_0)\right) ~\tau.$$ =

According to Bartczak & Breiter (2003), the potential of a significantly non-spherical ellipsoid can be efficiently approximated by means of a combination of two material segments: one of the length 2l₁ and mass $m_{\rm e} \mu_1$, placed symmetrically along the OX₁ axis, and one one with the imaginary length 2 l₃ = 2 i L₃ but with a real mass $m_{\rm e} \mu_3$, placed along the OX₃ axis. According to our recent experience, a more accurate and more convenient model can be obtained if a central point mass $m_{\rm e} \mu_0$ is added. Subtracting the Keplerian term from the potential of two segments and the central mass, one obtains

 

$$\mu \left[\frac{1-\mu_0}{R} - \frac{\mu_1}{2~l_1} ~ \ln{\left( \f... ...} \right)} - \frac{\mu_3}{L_3} ~\arctan\left( \frac{2 L_3}{s_3}\right) \right],$$ (40)

where $\mu = k^2 m_{\rm e} m_{\rm s}$,

  

$$\sqrt{R^2 + 2~X_1~l_1 + l_1^2}+ \sqrt{R^2 - 2~X_1~l_1 + l_1^2},$$                                        s₁ = (41)

$$\sqrt{2~\left( R^2 - L_3^2 + p_3\right)},$$ s₃ = (42)

$$\sqrt{(r^2- L_3^2)^2+4 X_3^2 L_3^2}.$$ p₃ = (43)

The gradient of V₁ is

  

$$\frac{\partial V_{1}}{\partial X_1} \mu~X_1 \left[~ \frac{2~\mu_1}{s_1~p_1} + \frac{2~\mu_3~s_3}{p_3~\left(s_3^2+4~L_3^2\right)} -\frac{1-\mu_0}{R^3} ~ \right],$$ = (44)

$$\frac{\partial V_{1}}{\partial X_2} \mu ~X_2 \left[ ~ \frac{2\mu_1~s_1}{p_1~\left(s_1^2-4~l_1^2\right... ...~\mu_3~s_3}{p_3~\left(s_3^2+4~L_3^2\right)} -\frac{1-\mu_0}{R^3} ~ \right],~~~~$$ = (45)

$$\frac{\partial V_{1}}{\partial X_3} \mu~X_3 \left[ ~ \frac{2\mu_1~s_1}{p_1~\left(s_1^2-4~l_1^2\right)} + \frac{2~\mu_3}{s_3~p_3} -\frac{1-\mu_0}{R^3} ~\right],$$ = (46)

where

$$p_1 = \sqrt{R^2 + 2~X_1~l_1 + l_1^2}~ \sqrt{R^2 - 2~X_1~l_1 + l_1^2}.$$ (47)

The parameters of our potential model can be easily determined by matching the coefficients of the Legendre series for V₁ with the C_0,0, C_2,0, C_2,2, C_4,0, and C_4,2 coefficients of the spherical harmonics expansion of a homogenous ellipsoid. A more detailed study of this model will be published in a separate article; here we only quote the formulas for the parameters

$$\sqrt{\frac{5~\left(a^2+b^2-2 c^2\right)}{7 }},$$                             l₁ = (48)

$$\sqrt{\frac{5~\left(a^2+3b^2-4 c^2\right)}{28 }},$$ L₃ = (49)

$$\mu_1 \frac{21}{25}~\frac{a^2-b^2}{a^2+b^2-2 c^2},$$ = (50)

$$\mu_3 \frac{84}{25}~\frac{b^2-c^2}{a^2+3 b^2-4 c^2},$$ = (51)

$$\mu_0 1- \mu_1 - \mu_3.$$ = (52)

Note that the real segment can be rejected ($\mu_1=0$) if the body has a=b, and the imaginary segment becomes negligible for b=c, when $\mu_3=0$.

Using the straight segments approximation is by no means a crucial feature. One may consider using the exact potential of an ellipsoid in terms of elliptic integrals or the spherical harmonics series, as long as the former is fast enough or the latter reasonably convergent.

   
4 Tangent maps

In many applications the knowledge of a trajectory $\vec{\zeta}(t)$ alone is not sufficient. The information how small variations of the initial conditions $\vec{\delta}(t)$ evolve in time is equally useful for the theoretical stability analysis as well as for practical problem of adjusting the modeled trajectory to observations. In principle, infinitesimal variations obey differential equations of motion and one might think about applying a splitting method to solve them. But the equations, albeit linear, are explicitly time-dependent due to the presence of the fiducial trajectory expressions $\vec{\zeta}(t)$ in their right-hand sides.

A practical solution to the problem of finding $\vec{\delta}(t)$ is based on the observation that symplectic (or Lie-Poisson) integration amounts to the iteration of a map $\Phi$, generating a sequence of values

$$\vec{\zeta}(t_0+n~h) = \vec{\zeta}_n = \Phi_h \vec{\zeta}_{n-1} = \Phi_h^n \vec{\zeta}_0.$$ (53)

If the initial conditions $\vec{\zeta}_0$ are perturbed by $\vec{\delta}_0$, a new sequence is generated, differing from $\vec{\zeta}_n$ by the evolving variations $\vec{\delta}_n$, namely

$$\vec{\zeta}_n + \vec{\delta}_n = \Phi_h^n (\vec{\zeta}_0 + \vec{\delta}_0).$$ (54)

The map $\Phi$ is usually nonlinear, so we have to approximate (54) using the linearized Taylor expansion - an approximation that becomes exact if $\vec{\delta}$ is considered infinitesimal. According to the chain rule for a differentiation of an n-fold composition of $\Phi$, we obtain

$$\Phi_h^n (\vec{\zeta}_0 + \vec{\delta}_0) = \Phi_h^n \vec{\... ... ~ {\bf A}_{n-1} \ldots {\bf A}_2 ~ {\bf A}_1 \vec{\delta}_0,$$ (55)

where ${\bf A}_j = (D \Phi_h)_j$ is the Jacobian matrix of $\Phi_h$evaluated at $\vec{\zeta}=\vec{\zeta}_{j-1}$. Substituting Eqs. (55) into (54), we obtain the evolution of $\vec{\delta}_n$ as the product of the Jacobian matrices

$$\vec{\delta}_n = {\bf A}_n \vec{\delta}_{n-1} = {\bf A}_n \ldots {\bf A}_1 \vec{\delta}_0.$$ (56)

Each matrix product is isomorphic with some linear transformation - the transformation that we call a tangent map. Thanks to the linearity of the tangent map, both sides of Eq. (56) can be multiplied by an arbitrary factor. It means that although formally $\vec{\delta}_j$ are considered either infinitesimal or "sufficiently small'', we may adopt any values of $\vec{\delta}_0$ in practical computation - a definite advantage with respect to the method of a direct integration of two nearby trajectories, where the roundoff error unavoidably follows the small displacements. Thus, we will rather speak about the tangent vector $\vec{\delta}$ instead of `variations'.

Let us now examine the expressions of tangent maps associated with $\Phi_h$ of our four-terms leapfrog. According to the definition (18), the tangent map matrix ${\bf A}_j$ is a product

 

$${\bf A}_j {\bf W}_{0,2,j} {\bf W}_{1,2,j}{\bf W}_{2,2,j}{\bf W}_{3,j} {\bf W}_{2,1,j} {\bf W}_{1,1,j} {\bf W}_{0,1,j},$$ (57)

where ${\bf W}_{i,k,j}$ is the Jacobian matrix of $\exp[\tau L_i]$ evaluated with the values of $\vec{\zeta}$that serve as the initial condition for the respective $\exp[\tau L_i]$. The second subscript has been added to indicate that the functional form of the matrices derived from the same hamiltonian $\mathcal{H}_i$ is similar, but ${\bf W}_{i,1,j}$ and ${\bf W}_{i,2,j}$ have different values due to the meanwhile evolution of $\vec{\zeta}$. In the following text we will skip the last two subscripts for the sake of brevity.

Another important remark concerns the actual form of tangent maps. Tangent maps are linear transformations and as such they can always be written as matrix products. But the matrix form is often too cumbersome and then it is easier to use the formulas resulting from the direct differentiation of the original nonlinear map.

4.1 Tangent map for $\mathcal{H}{_0}$

According to Sect. 3.1, $\exp[\tau L_0] = \Phi^{{\rm rot}}_\tau \circ \Phi^{{\rm Kep}}_\tau$, hence ${\bf W}_0 = {\bf W}^{\rm rot} {\bf W}^{\rm Kep}$. The tangent map for the Keplerian motion has been derived by Mikkola & Innanen (1999). Using their formulation in our case, when

$$\vec{\delta} = \left(% \begin{array}{c} \delta \vec{R} \\ \delta \vec{P} \\ \delta \vec{G} \\ \end{array}% \right),$$ (58)

requires a minor modification: all velocities $\vec{v}$ and their tangent vectors $\delta\vec{v}$should be replaced by $m^{-1}\vec{P}$ and $m^{-1}\delta\vec{P}$ respectively. The tangent vectors of angular momentum $\delta\vec{G}$ are constant in the Keplerian tangent map, because $\vec{G}$ itself is constant.

Straightforward differentiation of Eq. (29) leads to ${\bf W}^{\rm rot}$ - the tangent map for the rotation around $\vec{\hat{G}}$

$${\bf W}^{\rm rot} \vec{\delta} = \left(% \begin{array}{ccc... ...0} & {\bf0} & {\bf E} \\ \end{array}% \right) ~ \vec{\delta},$$ (59)

where

 

$$\delta {\bf M}_0 - \left( \sin{\theta} \left({\bf E} - \vec{\hat{G}}~\vec{\hat{G}}^{\rm T} \right) + \cos{\theta} ~{\bf Q}(\vec{\hat{G}}) \right) ~\delta \theta$$ =

$$+ 2 \sin^2{(\theta/2)} \left[ (\delta \vec{\hat{G}})~\vec{\hat{G}... ... \vec{\hat{G}})^{\rm T} \right] - \sin{\theta} ~ {\bf Q}(\delta \vec{\hat{G}}),$$ (60)

$$\delta \vec{\hat{G}} \left( \delta \vec{G} - \delta G ~ \vec{\hat{G}} \right)~ G^{-1},$$ = (61)

$$\delta G \left(\delta \vec{G} \right)^{\rm T} \vec{\hat{G}},$$ = (62)

$$\delta \theta \theta~(\delta G)~G^{-1}.$$ = (63)

In the dyadic notation of the above formulas, ${\bf B}=\vec{x} \vec{x}^{\rm T}$ stands for the square matrix with the elements B_ij = x_i x_j.

4.2 Tangent maps for $\mathcal{H}{_1}$ and $\mathcal{H}{_2}$

Tangent maps for the elementary rotations can be easily derived through the differentiation of the respective rotation matrices. Let

$${\bf M}'_1(\varphi_1) = \left(% \begin{array}{ccc} 0 & 0 &... ... -\cos{\varphi_1} & -\sin{\varphi_1} \\ \end{array}% \right),$$ (64)

and

$${\bf M}'_2(\varphi_2) = \left(% \begin{array}{ccc} -\sin{\... ...os{\varphi_2} & 0 & -\sin{\varphi_2} \\ \end{array}% \right).$$ (65)

Then, the tangent map due to $\mathcal{H}{_1}$ is

$${\bf W}_1 \vec{\delta} = \frac{\delta G_1~\varphi_1}{G_1} ... ... & {\bf0} & {\bf M}_1 \\ \end{array}% \right) ~ \vec{\delta},$$ (66)

and similarly

$${\bf W}_2 \vec{\delta} = \frac{\delta G_2~\varphi_2}{G_2} ... ... & {\bf0} & {\bf M}_2 \\ \end{array}% \right) ~ \vec{\delta}.$$ (67)

4.3 Tangent map for $\mathcal{H}{_3}$

The kick map (39) leads to the tangent map ${\bf W}_3$

$${\bf W}_3 \vec{\delta}= \left(% \begin{array}{ccc} {\bf E... ...t] & {\bf0} & {\bf E} \\ \end{array}% \right) ~ \vec{\delta},$$ (68)

where ${\bf V}_{RR}$ stands for the Hessian matrix composed of the second derivatives of V₁.

   
5 MLCE and MEGNO

Knowing the evolution of the tangent vector $\vec{\delta}$, we may study the qualitative properties of motion in the phase space by evaluating the Maximum Lyapunov Characteristic Exponent (MLCE) $\sigma$. For a continuous orbit $\vec{\zeta}(t)$ that satisfies the differential equation $\dot{\vec{\zeta}} = \vec{F}(\vec{\zeta})$, the MLCE is defined as

$$\sigma = \lim_{t \rightarrow \infty} \frac{1}{t} \ln \left... ...c{1}{t} \int_0^t \frac{\dot{\delta}(s)}{\delta(s)} {\rm d}s ,$$ (69)

where $\dot{\delta} = \delta^{-1} \vec{\delta}^{\rm T} (D\vec{F})^{\rm T} \vec{\delta}$ (Cincotta et al. 2003). For the iterations of a discrete time map $\vec{\zeta}_n = \Phi^n \vec{\zeta}_0$, with an associated tangent map $\vec{\delta}_n = (D\Phi)^n \vec{\delta}_0$, the definition of MLCE is

$$\sigma = \lim_{n \rightarrow \infty} \frac{1}{n} \sum_{k=1}^n \ln \left(\frac{\delta_k}{\delta_{k-1}}\right)\cdot$$ (70)

Roughly speaking, a positive value of MLCE implies a chaotic motion, which means that $\vec{\zeta}(t)$ and some other orbit $\vec{\zeta}'(t)$ that were infinitesimally close at t=0, diverge exponentially at the rate $\exp(\sigma t)$.

Evaluating the MLCE by means of numerical integration, one has to rely on some finite time estimate $\sigma(t)$ of $\sigma$. The convergence of $\sigma(t)$ is often slow and requires the integration on unreasonably long intervals. The Mean Exponential Growth factor of Nearby Orbits (MEGNO) function, invented by Cincotta & Simó (2000), is a good remedy for this difficulty. For a continuous orbit, the MEGNO function y(t)is defined as

$$y(t) = \frac{2}{t} ~ \int_{0}^{t} \frac{\dot{\delta}(s)}{\delta(s)} ~ s~ {\rm d}s.$$ (71)

The final conclusions about the nature of the studied orbit are drawn from the asymptotic behavior of the time-averaged MEGNO function Y(t)

$$Y(t) = \frac{1}{t} ~ \int_{0}^{t} y(s)~{\rm d}s.$$ (72)

For a discrete map, the MEGNO is computed from (Cincotta et al. 2003)

$$y(n) = \frac{2}{n} \sum_{k=1}^n k~\ln\left( \frac{\delta_k}{\delta_{k-1}}\right),$$ (73)

and

$$Y(n) = \frac{1}{n} \sum_{k=1}^n y(k).$$ (74)

Similarly to the MLCE case, we use a finite time estimate Y(t), but the behavior of MEGNO is quite different. Unlike $\sigma$, neither y nor Y have limits for a chaotic orbit. Actually, Y(t) asymptotically tends to

Y(t) = a t+b,

with a=0, b=2 for a quasi-periodic orbit, a=b=0 for a stable, isochronous periodic orbit, and $a = \frac{1}{2} \sigma,~b=0$for a chaotic orbit. Thus the MEGNO allows an indirect evaluation of MLCE in the chaotic case, but the convergence is much faster - usually it is enough to follow the interval of 10³ times the basic period of the system (unless a significant stickiness phenomenon appears).

The information that we obtain from the MEGNO in a non-chaotic case is more rich than what we get from the MLCE. From the formal standpoint, Y(t)converges to 2 for all kinds of quasi-periodic orbits; but the convergence rate is different depending on the local phase space environment and this provides an additional depth to the finite time estimators Y(t). Quasi-periodic orbits close to a stable periodic orbit tend to approach 2 from the left: the closer they are to the periodic orbit, the more slowly their Y(t) approach 2. Similarly, the neighborhood of an unstable periodic orbits slows down the convergence of Y(t) to 2 from the right. In these circumstances, the finite time (or even "short time'') estimators reveal the presence of periodic orbits between the grid points of the studied initial conditions or parameters.

A fixed step integrator, like the one we present in this paper, is equivalent to a discrete time map, so we have applied the definition of the MEGNO from Eqs. (73) and (74). The practical implementation of the algorithm was initially the following:

1.

After selecting the parameters and initial conditions $\vec{\zeta}_0$ and initializing y(0) = Y(0) = 0, the initial tangent vector $\vec{\delta}_0$ is generated with the components being random numbers from the (-0.5,0.5) range. The vector is then normalized in order to have $\delta_0=1$. This step was later modified as described below.

2.

Each time the nth step of the integrator is completed and the tangent vector $\vec{\delta}_n$ is computed, we update the values of y and Y according to

$$y(n) = \frac{n-1}{n}~y(n-1) + 2 ~ \ln \left(\frac{\delta_n}{\delta_{n-1}}\right),$$ (75)

and

$$Y(n) = \frac{(n-1) ~ Y(n-1) + y(n)}{n}\cdot$$ (76)

3.

After certain number of integrator steps (depending on the degree of chaoticity), the current tangent vector is renormalized to $\delta_n=1$ in order to avoid numerical overflow problems.

The random choice of $\vec{\delta}_0$ is important, because the convergence rate of Y(n) may significantly depend on the orientation of the initial tangent vector. If all initial conditions have the same initial tangent vector, than a double tadpole-shaped artifacts with small (even negative) MEGNO values appear, with the tadpoles' tails at the unstable equilibrium point and heads close to the stable equilibrium. According to our experience, this feature is generic whenever a libration zone appears on a phase plot (see Breiter et al. 2005, for a collection of such artifacts). The random choice of $\vec{\delta}_0$ usually leads to a slightly more "noisy'' plots of the MEGNO, but it removes such artificially created patterns - at least in principle.

But in the triaxial satellite limit of the Kinoshita problem we faced the difficulty, that the random choice of $\vec{\delta}_0$ components in the $(-\frac{1}{2},~\frac{1}{2})$ range did not lead to the removal of low MEGNO artifacts. Observing that ${\rm d}(\delta \vec{G})/{\rm d}t$is much smaller than ${\rm d}(\delta \vec{R})/{\rm d}t$ or ${\rm d}(\delta \vec{P})/{\rm d}t$, we realized, that even the random choice of $\vec{\delta}_0$ with all components equally treated actually preferred the direction of the slowest MEGNO growth. The problem was solved by the following device: at the beginning $\vec{\delta}_0^\ast$ is randomly generated with all components in the $(-\frac{1}{2},~\frac{1}{2})$ range, then one step of algorithm is performed propagating the values of $\vec{\delta}$; finally the lengths of the so obtained $\delta R^\ast$, $\delta P^\ast$, and $\delta G^\ast$are used as the weighting factors to generate actual $\vec{\delta}_0$with the components of $\delta\vec{R}_0$ equal to $\delta R^\ast$ times random numbers from the $(-\frac{1}{2},~\frac{1}{2})$ interval etc. This device practically removed systematic patterns of small Y(t).

   
6 Applications

6.1 The motion of Hyperion

After completing the formal tests of our program based on the algorithms described in previous sections, we chose the system of Saturn and its satellite Hyperion as the first realistic problem to be investigated. Wisdom et al. (1984) predicted the chaotic rotation of Hyperion by inspecting the phase space of the Kinoshita problem in the $(m_{\rm e}/m_{\rm s}) \rightarrow 0$limit. Their prediction was confirmed by the analysis of the photometric observations undertaken by Klavetter (1989a). Wisdom et al. (1984) studied mainly the planar case ($\vec{G}$ normal to the orbital plane and aligned with the axis of the maximum moment of inertia) that is perfectly suited for the analysis by means of the Poincaré section maps and reducible to an almost paradigmatic Beletsky equation (Beletsky 1965). But in the same paper, the authors notice and discuss the attitude instability effect arising when the planar restriction is suppressed. Due to this effect the widespread surface of section plot for Hyperion is actually misleading, because some regular regions of this plot prove to be chaotic in the spatial problem.

In order to study the phase space of Hyperion-like objects the following set of physical parameters was assumed (Klavetter 1989b): Hyperion's semi-axes $a=190~ \mbox{km}$, b=0.76316 a, c=0.6 a, its density $\rho = 1400 ~\mbox{kg}~{\rm m}^{-3}$, and Saturn's mass $m_{\rm s} = 5.685 \times 10^{26}~\mbox{kg}$. Initial conditions of Saturn in the Hyperion-fixed frame were derived from the Keplerian elements

$$a_{\rm s} = 1481.1 \times 10^3 ~ {\rm km},\quad e = 0.1, \quad \omega = M = 0,$$ (77)

and the initial momentum vector components were G₁ = G₂ = 0, and $G_3 = C~\nu$. In the definition of Keplerian elements we replaced velocities with momenta, which is the best way to handle the problem of two bodies in a rotating frame. We studied 144 761 orbits with various values of initial ascending node longitudes $\Omega$ from 0 to $180\hbox{$^\circ$ }$ with $0\hbox{$.\!\!^\circ$ }5$ grid, and rotation rates $\nu$ from 0 to 4 n with (n/100) grid, where $n = k~\sqrt{(m_{\rm s}+m_{\rm p})/a_{\rm s}^3}$stands for the mean motion of Saturn. The former of the variable parameters was $\Omega$, but due to $\omega=M=0$, its value determines the variable $\vartheta$ - the angle between the line of apsides and the longest axis of Hyperion.

For each orbit the MEGNO coefficient was evaluated on the interval of 10³ of orbital periods $2\pi/n$. The computations were performed twice and the two simulations differed only by the treatment of the tangent vectors. First, the initial variations $\vec{\delta}_0$ were selected with a planar restriction

$$\delta R_3 = \delta P_3 = \delta G_1 = \delta G_2 = 0.$$ (78)

We name this case RTM (Restricted Tangent Maps). Then, for the second run, all components of $\vec{\delta}_0$ were allowed to have nonzero values. We name this choice ATM (Arbitrary Tangent Map). Let it be emphasized, that both series contained the same orbits; the restrictions, distinguishing RTM from ATM, affected only the tangent vectors. The results are presented in Fig. 1. In all MEGNO plots the values of Y are classed into 9 bins labeled from A to I, as explained in Table 1.

In Fig. 1a we see the MEGNO map that exactly mimics the well know Hyperion surface of section pictures that can be found in papers like Wisdom et al. (1984), Klavetter (1989b), Black et al. (1995), or Breiter & Buciora (2000), not to mention numerous popular science texts. We recognize a chaotic sea with the main island of regular motion close to 1:1 spin-orbit resonance at $\vartheta = 0\hbox{$^\circ$ }$ and $\nu \approx 0.5~n$, surrounded by a small archipelago of a secondary resonance. The small island of 1:2 resonance is visible at $\vartheta = 90\hbox{$^\circ$ }$ and $\nu \approx 0.85~n$. Other larger islands surround the 9:4 resonance ( $\vartheta = 90\hbox{$^\circ$ }$, $\nu \approx 2~n$) and 2:1 spin-orbit coupling ( $\vartheta = 0\hbox{$^\circ$ }$, $\nu \approx 2.35~n$). The values of $\nu/n$ should not be confused with the resonant ratios p:q in question, because the vertical axes of all plots in this paper refer to instantaneous, initial values of angular rate that can differ from the mean angular rate by the current value of the sum of periodic terms. The best analogy at hand is the difference between the orbital angular rate at pericentre and the mean motion.

	Figure 1: MEGNO maps for Hyperion in RTM a) and ATM models b).
Open with DEXTER

According to Wisdom et al. (1984), the attitude instability brings chaos into the 1:1 and 1:2 resonance regions. Our MEGNO map in Fig. 1b confirms these results. Surprisingly, the 1:2 island is still visible, but now as a region of extremely high MEGNO (Y > 1000) with the Lyapunov times shorter than 10 days. The island has turned into an abyss. We also do not lose the trace of the 1:1 zone (including the small archipelago !) - this time as a "less chaotic'' region.

The stability of spin-orbit resonances depends on two factors: orbital eccentricity e and $\omega_0$ - a parameter that describes the departure of the primary's matrix of inertia from the prolate spheorid with A=B (Wisdom et al. 1984). For a homogeneous ellipsoid

$$\omega_0 = \sqrt{\frac{3~(B-A)}{C}} = \sqrt{3}~\sqrt{\frac{a^2-b^2}{a^2+b^2}},$$ (79)

and conversely

$$\frac{b}{a} = \sqrt{\frac{3-\omega_0^2}{3+\omega_0^2}}\cdot$$ (80)

Thus, regardless of the shortest semi-axis c, $\omega_0 =0$ when a=b. Using our values of semi-axes, we obtain for Hyperion $\omega_0 = 0.89$.

In order to demonstrate the dependence of resonant rotation states on $\omega_0$ within RTM and ATM models, we computed the MEGNO maps for 401 various rotation rates and 201 $\omega_0$ values spanning the range from 0 up to 1.2. The latter value is the limit obtained with b=c=0.6 a. Our simulations were performed for two initial values of $\vartheta$ equal 0 and $90\hbox{$^\circ$ }$. Thus the plots of Figs. 2a and 2c should be interpreted as a collection of the leftmost (or, due to symmetry, the rightmost) columns of plots similar to Fig. 1, but for different shapes of Hyperion. Analogously, Figs. 2b and 2d collect the middle columns of the plots similar to Fig. 1. Combining the information from $\vartheta = 0\hbox{$^\circ$ }$ and $\vartheta = 90\hbox{$^\circ$ }$ one may partially reconstruct the basic features of the whole phase space. The figures can be compared with Wisdom et al. (1984, Fig. 3).

Table 1: MEGNO map codes and Lyapunov time estimates for Hyperion.

	Figure 2: MEGNO for different shapes of Hyperion; RTM ( top) and ATM ( bottom).
Open with DEXTER

	Figure 3: MEGNO for different eccentricities of Hyperion orbit; RTM ( top) and ATM ( bottom).
Open with DEXTER

Figure 2 tells the story of bifurcations and resonance overlaps. The stable resonance curves are seen as strings of low MEGNO values. When they enter the chaotic sea, their width varies reflecting the width of a stable zone that surrounds the periodic orbit for a given value of $\omega_0$. From time to time the stable zone branches which is the evidence of a periodic orbit bifurcation. The branches terminate quickly because of the created secondary resonances overlap (Lichtenberg & Lieberman 1992). For a given p:q ratio, the resonance curve originates at $\nu/n = p/q$when $\omega_0 =0$, and extends first horizontally but then bending up or down as $\omega_0$ grows. The robustness of the 9:4 resonance in Figs. 2b and 2d is spectacular: it survives as long as the 1:2 case - up to $\omega_0 \approx 1.1$. Confronting the planar case (Figs. 2a,b) with the full problem (Figs. 2c,d), we notice the effects known from Fig. 1: low $\nu$ and 1:2 zones change their stable nature into an extreme chaos. There are also moderate chaos (type G) streaks in upper-right regions of Figs. 2c and 2d - apparently generated by some non-obvious resonance effect and definitely deserving future studies.

Figure 3 is constructed and arranged similarly to Fig. 2, but it shows the dependence of MEGNO at $\vartheta = 0\hbox{$^\circ$ }$ and $\vartheta = 90\hbox{$^\circ$ }$ as a function of orbital eccentricity of Hyperion (more precisely - of Saturn in our Hyperion-based reference frame). The value of $\omega_0$ was set to 0.89 according to the adopted ellipsoid of inertia. The maps at the top were computed with RTM; their counterparts at the bottom were obtained using nonrestricted tangent vectors. A sample of 201 eccentricity values was used, ranging from 0 to 0.4, allowing us to see the subsequent destruction of stable resonance areas after the chain of bifurcations. The influence of attitude instability is similar to previous cases. Interestingly, the 1:1 synchronicity becomes chaotic starting from very low eccentricities, and the 1:2 state is extremely chaotic even for the circular orbit.

	Figure 4: MEGNO maps for a Hyperion-like binary asteroid computed with RTM a) and ATM b) models.
Open with DEXTER

6.2 Hyperion-like binary asteroid

Is the attitude instability phenomenon intrinsically related to the negligible mass of the triaxial body? The simplest way to answer this question is to consider a hypothetic system consisting of Hyperion and a sphere with $m_{\rm s}$ of similar order of magnitude as $m_{\rm e}$. Decreasing the Saturn's mass we appropriately shrunk the orbital radius to conserve the the orbital mean motion. Thus we built a binary asteroid with a Hyperion-like ellipsoid and a sphere with the Hyperion's density and the radius of $250~{\rm km}$(i.e. $m_{\rm s}/m_{\rm e} \approx 5$). The orbital semi-axis $a_{\rm s}$was $800~{\rm km}$. Maintaining the mean motion, we can use directly the values of Lyapunov times associated with MEGNO bins in Table 1.

Mapping MEGNO as a function of initial $\vartheta$ and $\nu$, maintaining the eccentricity e=0.1, we obtained the results presented in Fig. 4 - the analogues of Fig. 1. Recognizing the characteristic 1:1 and 1:2 zones we notice that the 1:2 spin-orbit resonance reacts on passing from RTM to ATM similarly to the real Hyperion case and a similar chaotic wave is seen at low rotation rates. But, interestingly, the 1:1 synchronicity is not destabilized completely. Thanks to distinguishing different levels of chaoticity, we recognize a pattern resembling a strongly perturbed Second Fundamental Model of resonance Henrard & Lemaitre (1983).

	Figure 5: MEGNO maps for 90 Antiope in RTM a) and ATM models b).
Open with DEXTER

	Figure 6: MEGNO for different shapes of 90 Antiope; RTM ( top) and ATM ( bottom).
Open with DEXTER

6.3 Binary system of 90 Antiope

Thanks to the Keck Adaptive Optics observations (Merline et al. 2000) and extensive photometric data covering 5 oppositions (first 4 collected in Michaowski et al. 2004), the asteroid 90 Antiope is known to be an interesting example of a binary system with two comparable components. According to the recent results of Michaowski et al. (2004), Antiope is most likely an almost synchronous system with the rotation period of the primary $T_{\rm e} = 16\hbox{$.\!\!^{\rm h}$ }5047 \pm 0.0002$ and the orbital period of the satellite $T_{\rm s} = 16\hbox{$.\!\!^{\rm h}$ }5051 \pm 0.0002$. These results were obtained using a kinematic model with a spheroidal primary having semi-axes $a= 42.5~{\rm km}$, b=c=0.9 a and a spherical satellite with the radius of $42.5~{\rm km}$. A common density of the components was assumed as $\rho = 1300 ~{\rm kg}~{\rm m}^{-3}$ and the circular orbit of the satellite had the radius $a_{\rm s}= 170~{\rm km}$. In spite of the successful prediction of the observed lightcurves by means of the current model, its kinematic origin calls for a closer look from the dynamical point of view.

Figure 5 presents the dependence of MEGNO on the initial $\vartheta$ and $\nu$for Antiope. The values of MEGNO are coded according to Table 1, but to account for the difference of orbital periods, the estimates of Lapunov times from Table 1 should be divided by 30. The results obtained using RTM and ATM do not differ as significantly as in the previous examples. The 1:1 synchronicity zone is wide and stable in both cases, although the stochastic separatrix layer is quite thick. Few secondary resonances are also visible at $\nu/n \approx 1$ and $\vartheta$ between  $50\hbox{$^\circ$ }$ and  $70\hbox{$^\circ$ }$. Is it justified to claim that in a weakly spheroidal problem with comparable masses the attitude instability does not play any significant role? The answer is negative if we inspect Fig. 6. It is arranged similarly to Fig. 4: columns refer to the initial $\vartheta$ equal 0 (left) or $90\hbox{$^\circ$ }$ (right), and rows display RTM (top) or ATM (bottom). The appearance of secondary resonances bifurcating from the main 1:1 state is seen as the sequence of branchings in Fig. 6a. But comparing them with ATM map in Fig. 6c we see that three of the secondary resonances are attitude unstable - notably the one close to $\omega_0 \approx 0.37$. There is also the evidence of the high MEGNO wave at small $\omega_0$ and $\nu/n < 1$ that shrinks when the shape changes from a=b towards the b=c spheroid.

All previous examples discussed the planar Kinoshita problem. The last simulation scanned the whole range of orbital inclinations ( $0\hbox{$.\!\!^\circ$ }5$ grid) with respect to the equator of the spheroidal component, assuming circular orbital motion. The rotation rates were sampled in the $-3 \leqslant \nu/n \leqslant 3$ range with a step of 0.01. This time there was no reason to impose RTM, hence Fig. 7 consists of two parts only: 7a was obtained with the initial $\Omega = 270\hbox{$^\circ$ }$ and argument of latitude $M+\omega = 90\hbox{$^\circ$ }$, and 7b - with $\Omega = 0$ and $M+\omega = 90\hbox{$^\circ$ }$. The former leads to $\vartheta=0$ and the latter to $\vartheta = 90\hbox{$^\circ$ }$ when the inclination is 0.

The symmetry of both pictures with respect to their central points is obvious and it should be expected. The width of the 1:1 stable librations zone decreases with growing inclination (Fig. 7a) but at the same time the width of the 1:1 chaotic zone close to unstable equilibrium becomes smaller (Fig. 7b). As usually, the 1:2 resonance is the most capricious we can observe a sudden growth of the associated chaotic zone when inclination approaches $90\hbox{$^\circ$ }$. Of course, the same is true for 1:(-2) and retrograde orbits. Generally, there is a similarity between ($\nu/n,i$) and ( $-\nu/n,180\hbox{$^\circ$ }-i$) pairs. Polar orbits are chaotic for a slowly rotating primary if the orbital plane is perpendicular to the longest axis of the ellipsoid. If the orbital plane is aligned with the longest axis, polar orbits are stable. Interestingly, there are no vertical structures visible that might be attached to some critical inclination type phenomena.

	Figure 7: MEGNO as a function of initial inclination and rotation rate of 90 Antiope.
Open with DEXTER

7 Conclusions

Combining the MEGNO method with a fast integration algorithm results in a convenient tool for the study of spin-orbit interactions. The efficiency is gained thanks to few factors: the Lie-Poisson integrator allows to use a low-order method with relatively large stepsize without the immediate penalty of artificial energy drift and orthogonality loss, and MEGNO indicator has a faster convergence compared to traditional MLCE computation methods.

Although the applications presented in this paper were selected rather to illustrate the performance of the numerical algorithm, yet they reveal few interesting aspects of the studied problems. The difference between the values of MEGNO obtained for the same trajectory with restricted and arbitrary tangent maps is both understandable and intriguing. Suppose, that the trajectory is investigated by means of some Fourier spectrum method like the Frequency Analysis (Laskar 1993) or autocorrelation function; the conclusion should be similar to the RTM - no chaos. But the same trajectory investigated with Lyapunov exponents related methods without restrictions on tangent vectors qualifies the motion as chaotic. In our opinion, the discrepancy of results is due to the important difference between the chaos detectors based on the study of an orbit alone (Fourier type) an those referring to the divergence of nearby orbits (Lyapunov type).

Revisiting the motion of Hyperion we confirmed the results of Wisdom et al. (1984) concerning location of spin-orbit resonance states. We have also shown how the width of their stable zones varies depending on the moon's shape and orbital eccentricity.

For binary asteroids with comparable mass of components we have shown that chaotic tumbling is as common feature of the phase space as it is in the case of planetary satellites (Wisdom 1987). Any scenario of the evolution of binary asteroids should take this phenomenon into account regardless of the spin-up or spin-down effects included.

Our results concerning 90 Antiope contribute two important features to the current model of the system. First, the width of the regular libration zone surrounding the 1:1 synchronicity is too large to agree with the kinematic model of the tiny constant difference between the mean motion and rotation rate. On the other hand, the fact that no chaos has been observed in the lightcurve of Antiope may help to put better constraints on the shape of the primary, because some ratios of the ellipsoid's axes should be excluded according to Fig. 6c.

Acknowledgements

We thank Dr. Krzysztof Gozdziewski for numerous stimulating discussions about the theory and practice of the MEGNO. The research was carried under the Polish State Committee of Scientific Research (KBN) grant 1 P03D 020 27.

References

• Bartczak, P., & Breiter, S. 2003, Celest. Mech. Dynam. Astron., 86, 131 [NASA ADS] (In the text)
• Beletsky, V. V. 1965, Artificial Satellite Motion Relative to its Center of Mass, in Russian (Moscow: Nauka) (In the text)
• Black, G. J., Nicholson, P. D., & Thomas, P. C. 1995, Icarus, 117, 149 [NASA ADS] [CrossRef] (In the text)
• Breiter, S., & Buciora, M. 2000, Celest. Mech. Dynam. Astron., 77, 127 [NASA ADS] (In the text)
• Breiter, S., Wytrzyszczak, I., & Melendo, B. 2005, Adv. Space Res., in print (In the text)
• Cincotta, P. M., Giordano, C. M., & Simó, C. 2003, Physica D, 182, 151 [NASA ADS]
• Cincotta, P. M., & Simó, C. 2000, A&AS, 147, 205 [EDP Sciences] [NASA ADS] (In the text)
• Gozdziewski, K. 2003, A&A, 398, 1151 [EDP Sciences] [NASA ADS] [CrossRef] (In the text)
• Gozdziewski, K., Bois, E., Maciejewski, A. J., & Kiseleva-Eggleton, L. 2001, A&A, 378, 569 [EDP Sciences] [NASA ADS] [CrossRef] (In the text)
• Henrard, J., & Lemaitre, A. 1983, Celest. Mech., 30, 197 [NASA ADS] (In the text)
• Kinoshita, H. 1972, PASJ, 24, 423 [NASA ADS] (In the text)
• Klavetter, J. J. 1989a, AJ, 97, 570 [NASA ADS] [CrossRef] (In the text)
• Klavetter, J. J. 1989b, AJ, 98, 1855 [NASA ADS] [CrossRef] (In the text)
• Laskar, J. 1993, Celest. Mech. Dynam. Astron., 56, 191 [NASA ADS] (In the text)
• Lichtenberg, A. J., & Lieberman, M. A. 1992, Regular and Chaotic Dynamics, 2nd edn. (New York, Berlin, etc.: Springer), 258 (In the text)
• Maciejewski, A. J. 1996, Celest. Mech. Dynam. Astron., 63, 1 [NASA ADS]
• Marsden, J. E., & Ratiu, T. S. 1999, Texts in Applied Mathematics, Introduction to Mechanics and Symmetry (New York-Berlin-Heidelberg: Springer), 17, 291 (In the text)
• Merline, W. J., Close, L. M., Dumas, C., et al. 2000, BAAS, 32, 1309 (In the text)
• Merline, W. J., Weidenschilling, S. J., Durda, D. D., et al. 2002, Asteroids III, 289 (In the text)
• Michalowski, T., Bartczak, P., Velichko, F. P., et al. 2004, A&A, 423, 1159 [EDP Sciences] [NASA ADS] [CrossRef] (In the text)
• Mikkola, S., & Innanen, K. 1999, Celest. Mech. Dynam. Astron., 74, 59 [NASA ADS] (In the text)
• Olver, P. J. 1993, Graduate Texts in Mathematics, Application of Lie Groups to Differential Equations (New York-Berlin-Heidelberg: Springer), 107, 378
• Pavlov, A. I., & Maciejewski, A. J. 2003, Astron. Lett., 29, 552 [NASA ADS] [CrossRef] (In the text)
• Scheeres, D. J. 2002, Icarus, 159, 271 [NASA ADS] [CrossRef] (In the text)
• Touma, J., & Wisdom, J. 1994, AJ, 107, 1189 [NASA ADS] [CrossRef] (In the text)
• Wisdom, J. 1987, AJ, 94, 1350 [NASA ADS] [CrossRef] (In the text)
• Wisdom, J., & Holman, M. 1991, AJ, 102, 1528 [NASA ADS] [CrossRef] (In the text)
• Wisdom, J., Peale, S. J., & Mignard, F. 1984, Icarus, 58, 137 [NASA ADS] [CrossRef] (In the text)
• Yoshida, H. 1993, Celest. Mech. Dynam. Astron., 56, 27 [NASA ADS] (In the text)