EDP Sciences
Free access
Volume 506, Number 2, November I 2009
Page(s) 989 - 992
Section Celestial mechanics and astrometry
DOI http://dx.doi.org/10.1051/0004-6361/200811279
Published online 08 September 2009

A&A 506, 989-992 (2009)

On the predictability of unstable satellite motion around elongated celestial bodies
(Research Note)

E. Mysen

Institute of Theoretical Astrophysics, University of Oslo, PO Box 1029 Blindern, 0315 Oslo, Norway

Received 3 November 2008 / Accepted 28 August 2009

Context. Close satellite orbits around small and elongated celestial bodies can experience massive aperiodic changes in shape. According to a number of published works, the changes are continuous functions of the orbit elements of the satellite's unperturbed trajectory.
Aims. Later research, however, has revealed that the onset of instability is discrete, and highly correlated with the overlap of spin-orbit resonances. Since the interaction of resonances also can induce stochasticity in the orbiter's motion, we want to investigate more closely to what extent it is possible to predict the orbit evolution from the unperturbed elements.
Methods. Numerical simulations of a natural or artificial satellite's motion in a rotating gravity field of second order and degree are conducted using different algorithms and software.
Results. Consistent with the identification of resonance overlap as the responsible mechanism for the onset of orbit instability, we find that it is not always possible to predict qualitatively the outcome of a close encounter between satellite and the central body from the unperturbed orbit elements of the orbiter.
Conclusions. The massive aperiodic changes in orbit energy experienced by natural and artificial satellites in orbit around small elongated bodies exhibit properties characteristic for stochasticity.

Key words: space vehicules - celestial mechanics - methods: analytical - comets: individual: 67P/Churyumov-Gerasimenko

1 Introduction

Those who have simulated the motion of a satellite in orbit around a rotating and elongated central body, find that the orbit energy can experience massive changes in the vicinity of closest approach. From this observation, one can envision ejection of the orbiter from a bound trajectory, or collision with the central body.

According to Hamiltonian framework, spin-orbit resonances appear around the central body due to its rotation. Now, if we use the so-called resonance overlap criterion for chaos (Lichtenberg & Lieberman 1992), an aperiodic and unpredictable change in a satellite's orbit energy can be expected wherever the resonances overlap. The boundary between stable and unstable satellite orbits should, as a result, be almost discrete.

The overlap criterion has, however, been described as neither necessary nor sufficient for chaos (Lichtenberg & Lieberman 2002). Still, applications of the overlap criterion to this problem have yielded surprisingly precise predictions for the almost discrete transition between stable and unstable satellite orbits (Olsen 2006; Mysen et al. 2006; Mysen & Aksnes 2007). See also Petit et al. (1997).

It is therefore interesting that attempts to predict the changes in orbit energy analytically, by integration over a Keplerian trajectory, have been reported as successful, seemingly in contradiction with the identification of resonance overlap, and chaos, as the responsible mechanism for orbit instability. Also, the overlap criterion, with its almost discrete transition to instability, is not at terms with the continuous instability transition of the simple analytic approach.

The analytic technique has been applied a number of times to map the dynamical environment of asteroids and comets (Scheeres et al. 1996, 1998; Scheeres et al. 1999; Scheeres 2000a, 2003, 2004), and to predict changes in the rotation state of elongated objects during a planetary encounter (Scheeres 2000b; Scheeres 2001). Because of the wide spread applications of the approach, its importance to our idea of what can safely be done with a spacecraft in the vicinity of an elongated body, and its apparent inconsistencies with resonance overlap as responsible for orbit instability, we want to test the simple analytic approach more closely, especially since little numerical evidence has been given explicitly before.

In these tests of the analytic approch we have, in preparation for the ESA mission Rosetta, used parameters that correspond to close encounters between the spacecraft Rosetta and its primary target, comet 67P/Churyumov-Gerasimenko. However, it should be stressed that the considerations presented in this paper are general and therefore relevant to any close encounter of this kind.

2 Equations of motion

The details of any rigid body's gravity field can be expressed by a potential expanded in spherical harmonic functions (Heiskanen & Moritz 1967). To lowest order, the potential reads

V_1=\frac{\mu_{\rm c}}{r}~\left(\frac{r_{\rm c}}{r}\right)^2...
\end{displaymath} (1)

Here, $\mu_{\rm c}$ is the product of the constant of gravitation and the central body mass $m_{\rm c}=1\times{}10^{13}~{\rm kg}$, $r_{\rm c}\equiv{}2~{\rm km}$ is a reference radius, and r is the distance from the central body mass center to the satellite. The rotation of the central body is defined to be uniaxial, with period $P=12~{\rm h}$. The specific values are those of a nominal Rosetta target comet nucleus (see references in Mysen & Aksnes 2008). As for the constant gravity field coefficients, we use

\end{displaymath} (2)

which correspond to a fairly moderate central body elongation (Rossi et al. 1999). In Eq. (1), $P_{\rm 2m}(\sin\beta)$ are associated Legendre functions, where $\beta$ and $\lambda$ are the north latitude and east longitude, respectively, of the orbiter relative to the central body equator (body fixed xy-plane) and the so-called principal x-axis of the body (Goldstein 1980).

To obtain the corresponding acceleration $\ddot{\vec{r}}$ of the orbiter, one must calculate the gradient of the potential Eq. (1):

\end{displaymath} (3)

where Xi=(X,Y,Z) are reference system coordinates, illustrated together with some of the satellite's orbit elements in Fig. 1.
\end{figure} Figure 1:

Some angles of the satellite's orbit around the central body are shown. Above, i is the inclination, $\Omega $ the longitude of ascending node, and $u=\omega +f$ is the sum of the orbit's argument of pericenter, $\omega $, and the true anomaly, f.

Open with DEXTER
The spin axis of the central body is defined to be aligned with the reference Z-axis. For further details, see for instance Mysen & Aksnes (2005).

3 Predictability of unstable motion

It has previously been demonstrated that strong spin-orbit resonances exist in our problem, and that orbit shape instability is surprisingly well predicted by the resonance overlap criterion (see for instance Olsen 2006). Following the simplified first fundamental model of resonance, Mysen & Aksnes (2007) have parameterized the regions in satellite orbit semi-major axis, a, and eccentricity, e, space where overlap occurs. By definition, predictions based on the iterative integration over pure Keplerian motion fail in a resonance, and particularly so wherever the resonances overlap. Still, the structure of resonances and how they interact is highly non-trivial, and it is therefore possible that the previously mentioned simple analytic technique can be applicable.

3.1 Analytic approach

To explore this more closely, we look at the change in orbit energy from orbit apocenter, at time t0, to a time t1 when the semi-major axis of the orbit has stabilized (for instance the next apocenter):

\Delta{}F=\Delta{}\left(-\frac{\mu_{\rm c}}{2~a}+V_1\right)\...
...c{\mu_{\rm c}}{2}\left(\frac{1}{a}_1-\frac{1}{a_0}\right)\cdot
\end{displaymath} (4)

The change in energy can be calculated from basic Hamiltonian theory:

\Delta{}F=\int_{t_0}^{t_1}{\rm d}t~\frac{\partial{}F}{\partial{}t}\cdot
\end{displaymath} (5)

If the rotation is uniaxial, we have

...ht)^2~c_{22}~\frac{\partial{}Y^{\rm c}_{22}}{\partial{}\theta}
\end{displaymath} (6)

where $\theta$ is the rotational phase along the XY-plane, of the body fixed x-axis, as measured from the reference X-axis (Fig. 1). Note that $\dot{\theta}\equiv{}\omega_P=2\pi/P$, where P is the central body rotation period as before. Appropriate Fourier expansions of $Y^{\rm c}_{22}\equiv{}P_{22}(\sin\beta)\cos2\lambda$ can be found in Scheeres (2001) and Mysen et al. (2006), or in Kinoshita et al. (1974) for any degree and order of the gravity field. Accordingly

Y^{\rm c}_{22}=\sum_{q=-1}^{1}\Theta_q(i)~\cos{[~2(\theta-\Omega)-2q(f+\omega)~]}
\end{displaymath} (7)

with f as the true anomaly of the satellite in its orbit around the central body, and i as the orbit's inclination with respect to the central body equator, as before.

Already now, we will use a concept which is related the idea that we can say something about the orbit evolution from Keplerian motion: when $Y^{\rm c}_{22}$ above is substituted into Eq. (6), the q=-1 and q=0 components of Eq. (7) oscillate more rapidly than the q=1 term, so that the contribution from the former two to the aperiodic change in energy is negligible compared to the latter component. In addition, we do not assume that the inclination of the orbit with respect to the central body equator changes significantly during the change in orbit energy. These assumptions will be tested implicitly later on for some special cases. Finally, we have

$\displaystyle \Delta{}F_t\approx-2~\mu_{\rm c}~r_{\rm c}^2~c_{22}~\omega_P~\Theta_1(i)
\int_{t_0}^{t_1}{\rm d}t~\frac{\sin{2(\theta-\Omega{}-f-\omega)}}{r^3}$     (8)

where the sensitivity of the function

\end{displaymath} (9)

with respect to inclination defines what we mean by an insignificant change in i during pericenter passage.

The pioneer work of Scheeres (see references) then states that we can substitute ordinary time t in Eq. (8) with what we call pseudo-time t* that induces variation of the dynamical variables corresponding to Keplerian motion. Doing so results in

$\displaystyle \Delta{}F_{t^*}\approx-2~\mu_{\rm c}~r_{\rm c}^2~c_{22}~\omega_P~...
...left[~I_{\rm s}~\cos{2(\Omega+\omega)}-I_{\rm c}~\sin{2(\Omega+\omega)}~\right]$     (10)

$\displaystyle I_{\rm s}=\int_{t_0}^{t_1}{\rm d}t^*~\frac{\sin{2(\theta-f)}}{r^3}$     (11)
$\displaystyle I_{\rm c}=\int_{t_0}^{t_1}{\rm d}t^*~\frac{\cos{2(\theta-f)}}{r^3}\cdot$     (12)

These integrals will be evaluated numerically.

Under the assumption Eq. (10), the satellite orbit energy should not experience a significant aperiodic change if the longitude of pericenter, $\tilde{\omega}\equiv{}\Omega+\omega$, obeys

\tilde{\omega}_{{\rm min}}=\theta_p+k~\frac{\pi{}}{2}=\theta_0+\omega_P~\frac{T}{2}+k~\frac{\pi{}}{2}
\end{displaymath} (13)

where $k=\{...,-2,-1,0,1,2,..\}$. Above, T is the orbit period and $\theta_0$ and $\theta_{\rm p}$ are the rotational phases of the central body when the satellite is at orbit apocenter and pericenter, respectively. In the same way, we should expect a near maximum change in orbit energy when

\tilde{\omega}_{{\rm max}}(a)=\tilde{\omega}_{{\rm min}}(a)\pm{}\frac{\pi}{4}
\end{displaymath} (14)

if Eq. (10) holds.

3.2 Numerical simulations

We will now test if the truncation Eq. (8) can be an acceptable approximation of Eq. (4), and if the change in satellite orbit energy can be predicted using pseudo-time, inducing Keplerian motion, Eq. (10). For this purpose, the satellite orbit is initially defined to be almost directly prograde with respect to the central body's rotation, with initial semi-major axis $a(t_0)\equiv{}a_0=5~r_{\rm c}$. Here, it must be mentioned that outside the regions where the resonances overlap, the only quantity predicted by Eq. (10) is the difference in phase of the periodic variations in satellite orbit energy. That is, the approximation Eq. (10) does not ``clearly explain why particles in retrograde orbits about uniformly rotating body experience relatively small fluctuations in energy and angular momentum, as compared to particles in direct orbits'' (Scheeres 1999), see Eq. (9). Instead, it has been demonstrated (Olsen 2006; Mysen et al. 2006; Mysen & Aksnes 2007) that the stability of retrograde orbits is explained by a more fundamental principle than the amplitude of a slowly varying integrand being zero, namely that the resonances do not overlap.

In order to properly illustrate the strength of the interaction between satellite and central body in the vicinity of pericenter, we compare the corresponding changes in the orbit's semi-major axis a. The actual change $\Delta {}a$ is compared to $\Delta{}a_t$ derived from the truncated Eq. (8), and to the semi-major axis change $\Delta {}a_{t^*}$ from Eq. (10). These values are well approximated by

\begin{displaymath}\Delta{}a_{t,t^*}=\left(-\frac{2\Delta{}F_{t,t^*}}{\mu_{\rm c}}+\frac{1}{a_0}\right)^{-1}-a_0
\end{displaymath} (15)

where t0 represents the time when the satellite starts at first apocenter, and t1 is the time when the orbit energy stabilizes after closest approach. There is some arbitrariness to the choice of t1, but below, we have made sure that the conclusions are not affected by this choice. The results have been confirmed by two different softwares (Mathematica, Fortran 95) using different simulation algorithms (LSODA, Bulirsch-Stoer) and different precisions (15 and 32 digit mantissa).

Included as Fig. 2 are the semi-major axis changes for trajectories that are thought to maximize the change, Eq. (14), as functions of the initial satellite orbit's eccentricity e. Clearly, it can be seen that the truncated $\Delta{}a_t$, represented by one of the dashed curves, can serve as a good approximation of the actual change $\Delta {}a$. Also, the prediction based on the integration of the unperturbed Keplerian orbit is seen to work well for this maximum predicted from the analytic theory.

\end{figure} Figure 2:

The changes in semi-major axis as functions of the initial orbit's eccentricity e are shown. The actual change is given by $\Delta {}a$, while the prediction based on the integration over a Keplerian orbit is given by $\Delta {}a_{t^*}$. At first apocenter the semi-major axis is set to $a_0=5~r_{\rm c}$, and the orbit is defined to be almost directly prograde with respect to the central body rotation.

Open with DEXTER

Likewise, Fig. 3 shows how the analytic prediction compares to the actual change in semi-major axis when this change is predicted to be zero, Eq. (13). This time, there are discrepancies between the analytic prediction Eq. (10) and the actual change, which are clearly not caused by the truncation Eq. (4) since $\Delta{}a\approx{}\Delta{}a_{t}$.

\end{figure} Figure 3:

The actual change $\Delta {}a$ in semi-major axis of the satellite orbit during a close approach is compared to the analytic prediction $\Delta {}a_{t^*}$ for different eccentricities e of the initial orbit. The initial semi-major axis of the orbit is set to $a_0=5~r_{\rm c}$.

Open with DEXTER
Of course, it is not fair to compare the predictions of the simplified analytic technique to the actual change whenever the predicted change is close to zero, as such a comparison amplifies the roughness of the approach. However, a comparison with Fig. 2 demonstrates that the discrepancies are too large to be explained by the technique's roughness. Simulations with other initial semi-major axes indicate that the analytic approach performs well whenever the predicted change in orbit size is anticipated to be large, and not always so well whenever the change is predicted to be zero. The above result is a specification of the claim that the simple analytic technique agrees well with numerical simulations ``for high eccentricity elliptic orbits, parabolic orbits and hyperbolic orbits, but that accuracy begins to degrade sharply once the eccentricity of an elliptic orbit falls below a few tenths'' (Scheeres 1999).

3.3 Predictability

An important question remains: if the change in orbit energy is fairly well predicted by the analytic technique based on the integration over a Keplerian orbit, like in Fig. 2, is it possible then to predict the orbit evolution for successive close encounters? If so, it would apparently contradict the concept of resonance overlap as the mechanism responsible for the aperiodic change in orbit energy. One way to answer this question is to quantify the error $\delta{}a$ of the predicted change in semi-major axis during a close approach that alone, i.e. without regard to other satellite orbit parameters, would make the outcome of the next approach entirely uncertain. For instance, let the orbit's longitude of pericenter before first close encounter be $\tilde{\omega}_{{\rm min}}(a)$ of Eq. (13). The question is then: how much must the semi-major axis change during the encounter in order to turn a fixed $\tilde{\omega}$ from $\tilde{\omega}_{{\rm min}}(a)$ into $\tilde{\omega}_{{\rm max}}(a+\delta{}a)$, Eq. (14):

\begin{displaymath}0=\tilde{\omega}_{{\rm max}}(a+\delta{}a)-\tilde{\omega}_{{\r...
...t(\omega_P~\frac{\pi}{\sqrt{\mu_{\rm c}}}~a^{3/2}\right)-\pi/4


\begin{displaymath}\delta{}a\approx{}\frac{1}{6}~\frac{1}{\omega_P}~\sqrt{\frac{\mu_{\rm c}}{a_0}}\approx{}0.15~r_{\rm c}
\end{displaymath} (16)

where the specific value is for the adopted nominal Rosetta target and initial orbits of Figs. 2 and 3, namely $a_0=5~r_{\rm c}$. From Figs. 2 and 3, the error of the analytic technique is typically similar to or larger than $\delta{}a$ above. Hence, it is possible to consolidate the type of predictability presented by the analytic technique with the identification of resonance overlap as the underlying mechanism responsible for the rapid and aperiodic changes in orbit energy.

4 Discussion

Simulations show that the energy of a satellite orbit around an elongated central body can experience large aperiodic variations. Due to the apparent inconsistencies between an analytic technique used to predict these changes, and the previous identification of resonance overlap as the mechanism responsible for the observed effect, we have tested the analytic approach more closely.

It was found that the predictions provided by the studied analytic technique are limited to between zero and one closeencounter of the satellite with the central body. Accordingly, the predictions of the analytic technique were not found to be precise enough to rule out resonance overlap as the responsible mechanism for the observed effects.

The results strengthen the idea that spin-orbit resonances are essential when trying to understand the dynamical properties of and hazards for particles and spacecraft in orbit around asteroids and comets.

At last, it must be pointed out that we in the simulations sometimes see that the satellite is ejected from a bound orbit around the central body. The question is then if this behaviour is consistent with chaotic dynamics in its strictest sense since an orbit started in the chaotic regime is required to stay there. That is, an ejected satellite continues on its journey away from the central body, in a dynamical regime where two nearby trajectories do not diverge exponentially. Future studies will focus on this question.

The reason why the studied analytic technique works better for satellite trajectories which, according to the technique, should experience large changes in orbit energy during an encounter, is also of interest.

This work was possible due to support from the Research Council of Norway, project 184686.


  • Goldstein, H. 1980, Classical Mechanics (Reading: Addison-Wesley)
  • Heiskanen, W. A., & Moritz, H. 1967, Physical Geodesy (San Francisco: W.H. Freeman and Company)
  • Kinoshita, H., Hori, G., & Nakai, H. 1974, Modified Jacobi polynomials and its applications to expansions of disturbing functions, Annals of the Tokyo Astronomical Observatory, Second Series, Vol. XIV N.1
  • Lichtenberg, A. J., & Lieberman, M. A. 1992, Regular and Chaotic Dynamics (New York: Springer-Verlag)
  • Mysen, E., & Aksnes, K. 2005, A&A, 443, 691 [NASA ADS] [CrossRef] [EDP Sciences]
  • Mysen, E., & Aksnes, K. 2007, A&A, 470, 1193 [NASA ADS] [CrossRef] [EDP Sciences]
  • Mysen, E., & Aksnes, K. 2008, A&A, 490, 447 [NASA ADS] [CrossRef] [EDP Sciences]
  • Mysen, E., Olsen, Ø., & Aksnes, K. 2006, Planet. Space Sci., 54, 750 [NASA ADS] [CrossRef]
  • Olsen, Ø. 2006, A&A, 449, 821 [NASA ADS] [CrossRef] [EDP Sciences]
  • Petit, J.-M., Durda, D. D., Greenberg, R., Hurford, T. A., & Geissler, P. E. 1997, Icarus, 130, 1997 [CrossRef]
  • Rossi, A., Marzari, F., & Farinella, P. 1999, Earth Planets Space, 51, 1173 [NASA ADS]
  • Scheeres, D. J. 1999, Cel. Mech. Dyn. Astr., 73, 339 [NASA ADS] [CrossRef]
  • Scheeres, D. J. 2001, Cel. Mech. Dyn. Astr., 81, 39 [NASA ADS] [CrossRef]
  • Scheeres, D. J., Ostro, S. J., Hudson, R. S., & Werner, R. A. 1996, Icarus, 121, 67 [NASA ADS] [CrossRef]
  • Scheeres, D. J., Marzari, F., Tomasella, L., & Vanzani, V. 1998, Planet. Space Sci., 46, 649 [NASA ADS] [CrossRef]
  • Scheeres, D. J., Williams, B. G., & Miller, J. K. 2000a, J. Guid. Cont. Dyn., 23, 466 [CrossRef]
  • Scheeres, D. J., Ostro, S. J., Werner, R. A., Asphaug, E., & Hudson, R. S. 2000b, Icarus, 147, 106 [NASA ADS] [CrossRef]
  • Scheeres, D. J., Miller, J. K., & Yeomans, D. K. 2003, The Orbital Dynamics Environment of 433 Eros: A case study for Future Asteroid Missions, IPN Progress Report 42
  • Scheeres, D. J., Broschart, S., Ostro, S. J., & Benner, L. A. 2004, The Dynamical Environment About Asteroid 25143 Itokawa, ISTS 2004-d-36

All Figures

\end{figure} Figure 1:

Some angles of the satellite's orbit around the central body are shown. Above, i is the inclination, $\Omega $ the longitude of ascending node, and $u=\omega +f$ is the sum of the orbit's argument of pericenter, $\omega $, and the true anomaly, f.

Open with DEXTER
In the text

\end{figure} Figure 2:

The changes in semi-major axis as functions of the initial orbit's eccentricity e are shown. The actual change is given by $\Delta {}a$, while the prediction based on the integration over a Keplerian orbit is given by $\Delta {}a_{t^*}$. At first apocenter the semi-major axis is set to $a_0=5~r_{\rm c}$, and the orbit is defined to be almost directly prograde with respect to the central body rotation.

Open with DEXTER
In the text

\end{figure} Figure 3:

The actual change $\Delta {}a$ in semi-major axis of the satellite orbit during a close approach is compared to the analytic prediction $\Delta {}a_{t^*}$ for different eccentricities e of the initial orbit. The initial semi-major axis of the orbit is set to $a_0=5~r_{\rm c}$.

Open with DEXTER
In the text

Copyright ESO 2009