A&A 409, 1135-1140 (2003)
DOI: 10.1051/0004-6361:20031162

Canonical modelling of coorbital motion in Hill's problem using epicyclic orbital elements

P. Gurfil - N. J. Kasdin

Mechanical and Aerospace Engineering Department, Princeton University, Princeton NJ 08544, USA

Received 3 June 2003 / Accepted 23 June 2003

Abstract
This paper presents a Hamiltonian approach to modelling relative coorbital motion based on derivation of canonical coordinates for the relative coorbital dynamics. The Hamiltonian formulation facilitates the modelling of high-order terms and orbital perturbations while allowing us to obtain closed-form solutions to the relative coorbital motion in Hill's restricted three-body problem. First, the Hamiltonian is partitioned into a linear term and a high-order term. The Hamilton-Jacobi equations are solved for the linear part by separation, and new constants for the relative motions are obtained, called epicyclic orbital elements. The influence of the gravitational interaction between the coorbiting satellites is incorporated into the analysis by a variation of parameters procedure.

Key words: celestial mechanics - planets and satellites: general - reference systems

1 Introduction

The motion of any number of mutually gravitating satellites, confined to follow the same mean orbit about a massive primary, is usually referred to as coorbital motion. The coorbital motion problem can be modelled by using either a restricted-three body approach, in which one of the coorbiting objects is assumed to have an infinitesimal mass (Libre & Ollè 2001), or a more general treatment in which the coorbital satellies are assumed to have small, but nonzero, mass (Hénon 1969; Hénon & Petit 1986). The latter approach was originally developed by Hill (1878) and is known as Hill's restricted three-body problem. Applications of Hill's problem are ubiquitous (see e.g. Vilac & Scheers 2003; Scheeres et al. 2003, and the references therein). It has been shown that Hill's problem is characterized by the same generality as that of the restricted three-body problem, in the sense that the mass ration of the coorbiting objects is arbitrary (Hénon 1986).

The restricted three-body problem yields a very rich dynamical structure. In this work, we shall be primarily interested, however, in the 1:1 mean motion commensurability. For this problem, there are a few topologically distinct orbits: Tadpole orbits, taking place about the triangular Lagrangian equilibrium points L₄ and L₅ (e.g. Trojan asteroids); Horseshoe orbits about the Lagrangian points L₃-L₄-L₅ (e.g. the Saturnian satellites Epimetheus and Janus); and retrograde/prograde coorbital motion (e.g. rings of Saturn).

A number of authors have addressed the coorbital motion problem. Libre & Ollé (2001) have utilized the restricted three-body problem framework to characterize the coorbital motion of Saturnian satellites. Yoder et al. (1983) have offered a simplified qualitative framework for modelling coorbital motion. In a later work, Salo & Yoder (1988) gave sufficiency conditions for the stability of a system of Ncoorbital objects. Namouni (1999) performed an extensive study of coorbital motion using Hill's three-body problem and an orbital elements-based generating solution, an approach reminiscent to that of Winter & Murray (1997), who utilized orbital elements coupled with Lagrange's planetary equations in order to study resonant motion. Broucke (1999) derived a Lagrangian formulation for the study of the 1:1 motion commensurability in the restricted three-body problem. However, thus far, the literature has not presented a canonical modelling methodology for the coorbital relative motion problem.

This paper is aimed at developing canonical orbital elements for modelling the relative coorbital motion problem using Hill's equations. In other words, we attempt to find Delaunay-like canonical elements for the relative motion.

In order to solve the problem, we first write the Lagrangian for the coorbital motion, and then perform a Legendre transformation to find the Hamiltonian. We solve the Hamilton-Jacobi equations and treat the gravitational interaction of the coorbiting objects as a perturbation. The new orbital elements, which we termed epicyclic orbital elements, are constants of the relative motion. Due to the fact that the epicyclic orbital elements are canonical, any given perturbation can be modelled using Hamilton's equation. This methodology offers a simple and general framework for modelling coorbital motion. We illustrate the use of the newly developed methodology by characterizing a few periodic orbits using a numerical search procedure.

2 The Lagrangian for coorbital motion

We consider the coorbital motion of two bodies of masses m₁ and m₂ about a primary of mass M, with $m_1,m_2\ll M$. The coorbital motion can be described by any number of possible coordinate systems. The best coordinate system for our purposes is the one in which the Hamilton-Jacobi equation most easily separates. To this end, Cartesian coordinates turn out to be the most convenient. The work in this paper will be thus confined to a rotating Cartesian Hill frame. This coordinate system, denoted by $\mathscr{R}$, is defined by the unit vectors $\hat{x},~\hat{y},~\hat{z}$. The origin of this coordinate system is set on a circular reference orbit of radius a about the primary. It is rotating with mean motion $n=\sqrt{GM/a^3}$. The reference orbit plane is the fundamental plane, the positive $\hat{x}$-axis points radially outward, the $\hat{y}$-axis is rotated $90 \deg$ in the direction of motion and lies in the fundamental plane, and the $\hat{z}$-axis completes the setup to yield a Cartesian dextral system.

For simplicity, we first treat the case of coorbital motion with respect to a circular reference orbit. This is the most common problem and should easily reduce to Hill's equations. We start with this case because of its simplicity, allowing us to focus attention on the details of the method. Nevertheless, we find that the resulting canonical perturbation equations still provide new and meaningful results. In future work we will present the more involved case of arbitrary elliptical orbits.

We initially assume zero gravitational interaction between the coorbital bodies. We then introduce the gravitational interaction as generalized forces in the Lagrangian formulation and an interaction potential in the Hamiltonian formulation.

The first step is to develop the Lagrangian of the coorbital relative motion in the rotating frame $\mathscr{R}$. The velocities of the coorbiting objects in $\mathscr{R}$ is given by the usual equation:

$$\vec{v}_i={^\mathscr{I}\mathbf{\mbox{\boldmath$\omega$ }}_i^... ...thscr{I}\bfomega_i^\mathscr{R}}\times\mbox{\boldmath$\rho$ }_i$$ (1)

where $\vec r\in\mathbb R^3$ is the inertial position vector along the reference orbit, $\mbox{\boldmath$\rho$ }_i=[x_i,~y_i,~z_i]^T\in\mathbb R^3$is the relative position vector of each object in the rotating frame, and ${^\mathscr{I}\mbox{\boldmath$\omega$ }^\mathscr{R}}=[0,~0,~n]^T$ is the angular velocity of the rotating frame $\mathscr{R}$ with respect to an inertial frame $\mathscr{I}$. Assuming a circular reference orbit, denoting $\Vert\vec {r}\Vert=a$ and substituting into Eq. (1), we can write the velocity of each coorbiting object in a component-wise notation:

$$\vec{v}_i=\left [ \begin{array}{c} \dot{x}_i-ny_i \\ \dot{y}_i+nx_i+na \\ \dot{z}_i \end{array} \right ].$$ (2)

The kinetic energy per unit mass is given by

$$\mathcal K_i = \frac{1}{2}\Vert\vec v_i\Vert^2.$$ (3)

The potential energy (due to a spherical primary) per unit mass of the coorbiting objects, whose position vectors are denoted by $\vec{R}_i$, is the usual specific gravitational potential written in terms of $\rho_i=\Vert\mbox{\boldmath$\rho$ }_i\Vert$ and expanded using Legendre polynomials:

$$\mathcal U_i -\frac{GM}{\Vert\vec{R}_i\Vert}=-\frac{GM}{\Vert \vec{r}+\mbox{\b... ...\left ( \frac{\displaystyle \rho_i}{\displaystyle a} \right )^2 \right ]^{1/2}}$$ =

$$-\frac{GM}{a} \sum_{k=0}^{\infty} P_k(\cos\alpha_i)\left(\frac{\rho_i}{a} \right)^k$$ = (4)

where the $P_k(\cos\alpha_i)$ are the Legendre polynomials,

$$\cos\alpha_i = -\frac{\mbox{\boldmath$\rho$ }_i\cdot\vec{r}}{a\rho_i} = \frac{-x_i}{\sqrt{x_i^2+y_i^2+z_i^2}}$$ (5)

and $\alpha_i$ is the angle between the reference orbit radius vector and the position vector of each coorbiting object relative to the reference orbit.

The Lagrangian $\mathcal L$ is now easily found by subtracting the total potential energy from the total kinetic energy,

$$\mathcal{L}=\frac{1}{2}\sum_{i=1}^2 \left \{ (\dot{x}_i-ny_i)... ...nfty} P_k(\cos\alpha_i)\left(\frac{\rho_i}{a} \right )^k \cdot$$ (6)

As in the treatment leading to Hill's equations for relative motion, we examine only small deviations from the reference orbit. Thus, we only consider the first three terms of the potential energy for each orbiting body,

$$\mathcal U^{(0)}_i=-\frac{GM}{a}-\frac{GM}{a^2}\rho_i\cos\alp... ...}\rho_i^2\left (\frac{3}{2}\cos^2\alpha_i-\frac{1}{2} \right )$$ (7)

and then use Eq. (5) to find the low order Lagrangian,

$$\mathcal{L}^{(0)}=\sum_{i=1}^2\frac{1}{2}\left (\dot{x_i}^2+\... ...ht ) +\frac{3}{2}n^2a^2+\frac{3}{2}n^2x_i^2-\frac{n^2}{2}z_i^2$$ (8)

with the perturbed part of the Lagrangian equal to the higher order terms in the potential [ $\mathcal{O}((\rho_i/a)^3)$]. Hill's equations for the relative motion of the two satellites in the rotating frame are obtained by subtracting the Euler-Lagrange equations (with the Lagrangian given by Eq. (8)) for each orbiting body using the generalized coordinates $\vec{q}_i=\mbox{\boldmath$\rho$ }_i$,

$$\frac{\rm d}{{\rm d}t}\left(\frac{\partial\mathcal L^{(0)}... ...\frac{\partial\mathcal L^{(0)}}{\partial\vec{q}_2}\right ]=0,$$ (9)

resulting in Hill's equations,

   

$$\ddot{x}-2n\dot{y}-3n^2x = -\frac{\partial\Psi}{\partial x},$$ (10)

$$\ddot{y}+2n\dot{x} = -\frac{\partial\Psi}{\partial y},$$ (11)

$$\ddot{z}+n^2z = -\frac{\partial\Psi}{\partial z},$$ (12)

where x = x₁-x₂, y = y₁-y₂, z = z₁-z₂ and $\Psi$ is the interaction potential, given by

$$\Psi = -\frac{G(m_1+m_2)}{r}=-\frac{n^2 a^3\mu}{r},$$ (13)

with $r=\sqrt{x^2+y^2+z^2}$, and

$$\mu=\frac{m_1+m_2}{M}\cdot$$ (14)

A similar result to that of Eqs. (10)-(12) is obtained by expanding the nonlinear equations of motion with respect to the mass parameter (Hill 1878):

$$\epsilon = \left (\frac{\mu}{3}\right)^{\frac{1}{3}}\cdot$$ (15)

The equilibria of Eqs. (10)-(12) are located at $(x=\epsilon,~y=0,~z=0)$ and $(x=-\epsilon,~y=0,~z=0)$, and correspond to the Lagrangian points L₂ and L₁, respectively.

Based on Eqs. (10)-(12), using the new generalized coordinate $\vec{q}=\vec{q}_1-\vec{q}_2=(x,~y,~z)$, we can define the relative motion Lagrangian

$$\mathcal{L}^{(0)}_{\rm r}=\frac{1}{2}\left (\dot{x}^2+\dot{y}... ...ht ) +\frac{3}{2}n^2a^2 +\frac{3}{2}n^2x^2-\frac{n^2}{2}z^2,$$ (16)

and proceed with the analysis using this Lagrangian.

3 The Hamilton-Jacobi solution

Our overall objective is to find the Hamiltonian of the relative motion and divide it into a linear part and a perturbed part,

$$\mathcal{H}_{\rm r}=\mathcal{H}^{(0)}_{\rm r} + \mathcal{H}^{(1)}_{\rm r}$$

and then solve the Hamilton-Jacobi equation for the linear system without the interaction potential. This solution will provide us with new canonical coordinates and momenta that are constants of the relative motion. The perturbation equations will then show how these constants vary under the gravitational interaction between the coorbiting objects.

Finding the Hamiltonian for the system is straightforward. First, we utilize the normalization n=a=1. The canonical momenta are then found from the usual definition:

 

$$\frac{\partial\mathcal{L}^{(0)}_{\rm r}}{\partial\dot{x}} = \dot{x}-y$$                   p_x =

$$\frac{\partial\mathcal{L}^{(0)}_{\rm r}}{\partial\dot{y}} = \dot{y}+x+1$$ p_y = (17)

$$\frac{\partial\mathcal{L}^{(0)}_{\rm r}}{\partial\dot{z}} = \dot{z}$$ p_z =

and then, using the Legendre transformation, the unperturbed Hamiltonian for relative motion in Cartesian coordinates is found:

$$\mathcal{H}^{(0)}_{\rm r}=\frac{1}{2}(p_x+y)^2+\frac{1}{2}(p_... ...2 +\frac{1}{2}p_z^2-\frac{3}{2}-\frac{3}{2}x^2+\frac{1}{2}z^2.$$ (18)

This Hamiltonian is used to solve the Hamilton-Jacobi equation (Goldstein 1980). Because the Hamiltonian is a constant, Hamilton's principal function easily separates into a time dependent part summed with Hamilton's characteristic function,

$$S(x,y,z,t)=W(x,y,z)-\alpha'_1 t$$

where $\alpha'_1$ is the constant value of the unperturbed Hamiltonian, $\mathcal{H}^{(0)}_{\rm r}$. The Hamilton-Jacobi equation then reduces to a partial differential equation in W(x,y,z):

$$\frac{1}{2}\left (\frac{\partial W}{\partial x}+y\right )^2+\... ...ight )^2-\frac{3}{2}-\frac{3}{2}x^2+\frac{1}{2}z^2=\alpha'_1 .$$ (19)

Not unexpectedly, the z-coordinate easily separates. Separating the characteristic function as W(x,y,z)=W'(x,y)+W₃(z), the HJ equation separates into:

  

$$\frac{1}{2}\left (\frac{{\rm d}W_3}{{\rm d}z} \right )^2 + \frac{1}{2}z^2 = \alpha_2$$ (20)

$$\frac{1}{2}\left (\frac{\partial W'}{\partial x}+y\right )^2+\fra... ... W'}{\partial y}-x-1\right )^2-\frac{3}{2}x^2 = \alpha'_1+\frac{3}{2}-\alpha_2.$$ (21)

Equation (20) is just the HJ equation for simple harmonic motion. It is easily solved via quadrature:

 

$$\int \sqrt{2\alpha_2-z^2}{\rm d}z$$                       W₃(z) =

$$\frac{1}{2} \left [ z\sqrt{2\alpha_2-z^2}+2\alpha_2\sin^{-1}\left(\frac{z}{\sqrt{2\alpha_2}} \right ) \right ]\cdot$$ = (22)

The solution of Eq. (21) is more subtle. We separate by using the well known constant of integration of the CW equations. Setting $\alpha_3$ equal to the integral of Eq. (11) with $\Psi=0$ and putting it in terms of the canonical momentum, the third integration constant of the HJ equation is:

$$\alpha_3=p_y+x-1.$$ (23)

Using the fact that $p_y=\partial W'/\partial y$, the remaining HJ equation (Eq. (21)) separates if we let

W'(x,y)=W₁(x)+W₂(y)-yx (24)

so that

$$\frac{{\rm d}W_2}{{\rm d}y} = \alpha_3 +1$$

and thus $W_2 = (\alpha_3+1)y$ by quadrature. Equation (21) then simplifies to:

$$\left (\frac{{\rm d}W_1}{{\rm d}x}\right )^2+x^2-4\alpha_3x=2\alpha_1-\alpha_3^2$$ (25)

where we have used $\alpha_1=\alpha_1'-\alpha_2+3/2$. This equation is again easily integrated for W₁ by quadrature,

$$W_1=\int \sqrt{2\alpha_1-\alpha_3^2+4\alpha_3x-x^2}{\rm d}x$$ (26)

yielding

$$W_1=\left(\frac{x}{2}-\alpha_3\right)\sqrt{2\alpha1-\alpha_3^... ...lpha_3}{\sqrt{2\alpha1-\alpha_3^2+4\alpha_3x-x^2}}\right)\cdot$$ (27)

The final generating function from the solution of the low-order HJ equation is thus given by:

$$S(x,y,z,\alpha_1,\alpha_2,\alpha_3,t)=W_1(x)+W_2(y)+W_3(z) -(\alpha_1 +\alpha_2) t.$$ (28)

(Note that we have omitted the constant 3/2 as it does not affect the solution). It is straightforward to express the new canonical momenta $(\alpha_1,\alpha_2,\alpha_3)$ in terms of the original Cartesian positions and velocities (and thus in terms of the initial conditions). For instance, $\alpha_3$ is given by Eq. (23) using Eq. (17). Equation (20) is used to find $\alpha_2$, substituting p_zfrom Eq. (17) for ${\rm d}W_3/{\rm d}z$. Finally, $\alpha_1=\alpha'_1-\alpha_2+3/2$ is simply the value of the Hamiltonian and is thus given by Eq. (18) with the momenta substituted from Eq. (17). The result is:

   

$$\alpha_1 \frac{1}{2}\dot{x}^2 +\frac{1}{2}\dot{y}^2-\frac{3}{2}x^2$$ = (29)

$$\alpha_2 \frac{1}{2}\dot{z}^2+\frac{1}{2}z^2$$ = (30)

$$\alpha_3 \dot{y} + 2 x.$$ = (31)

The canonical coordinates ( $\beta_1,\beta_2,\beta_3$) are found via the partial derivatives of the generating functions in Eq. (28) with respect to each of the new canonical momenta,

$$\beta_i=\frac{\partial \left[S(x,y,z,\alpha_1,\alpha_2,\alpha_3,t)+\alpha'_1t\right]}{\partial\alpha_i}$$ (32)

yielding

   

$$\beta_1 -\tan^{-1}\left(\frac{3x+2\dot{y}}{\vert\dot{x}\vert}\right)$$ = (33)

$$\beta_2 \tan^{-1}\left(\frac{z}{\vert\dot{z}\vert}\right)$$ = (34)

$$\beta_3 -(3\dot{y}+6x)\tan^{-1}\left(\frac{3x+2\dot{y}}{\vert\dot{x}\vert}\right)-2\dot{x}+y.$$ = (35)

Solving Eqs. (29)-(35) for $x,~y,~\textrm{and}\:z$ yields the generating solution for the Cartesian relative position components in terms of the new constants of the motion, the canonical momenta ( $\alpha_1, \alpha_2, \alpha_3$) and the canonical coordinates ( $\beta_1,\beta_2,\beta_3$):

   

$$2\alpha_3 + \sqrt{2\alpha_1+3\alpha_3^2}\sin(t+\beta_1)$$                                        x(t) = (36)

$$\beta_3 - 3\alpha_3 (t+\beta_1) + 2\sqrt{2\alpha_1 + 3\alpha_3^2} \cos(t+\beta_1)$$ y(t) = (37)

$$\sqrt{2\alpha_2} \sin(t+\beta_2).$$ z(t) = (38)

From Eqs. (29)-(35) (or, alternatively, by differentiating Eqs. (36)-(37) with respect to time) we can also obtain the expressions for the Cartesian relative velocity components in terms of $\alpha_1, \alpha_2, \alpha_3$ and $\beta_1,\beta_2,\beta_3$:

   

$$\dot{x}(t) \sqrt{2\alpha_1+3\alpha_3^2}\cos(t+\beta_1)$$ = (39)

$$\dot{y}(t) - 3\alpha_3 - 2\sqrt{2\alpha_1 + 3\alpha_3^2} \sin(t+\beta_1)$$ = (40)

$$\dot{z}(t) \sqrt{2\alpha_2} \cos(t+\beta_2).$$ = (41)

Thus, as is well known from solutions of Hill's unperturbed equations, the motion consists of a periodic out-of-plane oscillation parameterized by $\alpha_2, \beta_2$, a periodic in-plane motion described by $\alpha_1, \beta_1$, and a secular drift in y given by $\alpha_3$. The well known invariance with y is given by the shift $\beta_3$.

We call the new constants of the motion $\Xi=[\alpha_1,~\alpha_2,~\alpha_3,~\beta_1,~\beta_2,~\beta_3]$epicyclic orbital elements for the relative motion. They are defined on the manifold $\mathscr{O}\times\mathbb S^3$, where $\mathscr{O}=\mathbb R\times\mathbb R_{\ge 0}\times\mathbb R\subset\mathbb R^3$.

4 Modified epicyclic elements

The epicyclic elements described above provide a convenient parametrization of a first-order relative motion orbit in terms of amplitude and phase. However, variational equations presented later for these elements to account for the gravitational interaction can become quite complicated and numerically sensitive. This is particularly a concern when some of the amplitudes approach zero, resulting in the phase terms becoming ill-defined. For these situations it is convenient to introduce an alternative set of constants in terms of amplitude variables only. We call these modified epicyclic orbital elements, label them $\Xi' = \left [\alpha'_1,\alpha'_2,\alpha'_3,\beta'_1,\beta'_2,\beta'_3\right]$, and find them via the point transformation:

      

$$\alpha'_1 \sqrt{2\alpha_1+3\alpha_3^2}\cos\beta_1$$ = (42)

$$\beta'_1 \sqrt{2\alpha_1+3\alpha_3^2}\sin\beta_1$$ = (43)

$$\alpha'_2 \sqrt{2\alpha_2}\cos\beta_2$$ = (44)

$$\beta'_2 \sqrt{2\alpha_2}\sin\beta_2$$ = (45)

$$\alpha'_3 \alpha_3$$ = (46)

$$\beta'_3 \beta_3 - 3\alpha_3\beta_1.$$ = (47)

It can easily be shown that the transformation in Eqs. (42)-(47) is symplectic. That is, letting $M=\partial\Xi'/\partial\Xi$, it is straightforward to show that:

MJM^T = J (48)

where J is the symplectic matrix [0,I;-I,0]. Thus, the new elements are also canonical and satisfy Hamilton's equations. In some cases, the variational equations for these elements will be much easier to work with. In terms of the new canonical momenta $(\alpha'_1,\alpha'_2,\alpha'_3)$ and new canonical coordinates $(\beta'_1,\beta'_2,\beta'_3)$, the generating solution becomes:

   

$$2\alpha'_3 + \alpha'_1\sin(t) + \beta'_1\cos(t)$$                         x(t) = (49)

$$\beta'_3 - 3\alpha'_3 t - 2\beta'_1\sin(t) + 2\alpha'_1 \cos(t)$$ y(t) = (50)

$$\beta'_2\cos(t) + \alpha'_2\sin(t).$$ z(t) = (51)

Variations of the parameters $\Xi$ or $\Xi'$ due to perturbations, such as gravitational interaction, can be obtained via Hamilton's equations on a perturbing Hamiltonian, and the resulting time varying parameters $\Xi(t)$ or $\Xi'(t)$ can then be substituted into the generating solutions (36)-(38) or (49)-(51) to yield the exact relative motion description in the configuration space $\mathbb{R}^3$. This is the subject of the next section.

5 Gravitational interaction analysis via variation of parameters

The primary value of the canonical approach is the ease with which equations for the variations of parameters can be found. For example, the variations of the epicyclic orbital elements are given by Hamilton's equations on the perturbation Hamiltonian, $\mathcal{H}^{(1)}_{\rm r}$:

  

$$\dot{\beta}'_i \frac{\partial \mathcal{H}^{(1)}_{\rm r}}{\partial \alpha'_i}=\fr... ...tial\mathcal{H}^{(1)}_{\rm r}}{\partial z}\frac{\partial z}{\partial \alpha'_i}$$ (52)

$$\dot{\alpha}'_i -\frac{\partial \mathcal{H}^{(1)}_{\rm r}}{\partial \beta'_i}=-\l... ...{H}^{(1)}_{\rm r}}{\partial z}\frac{\partial z}{\partial \beta'_i}\right) \cdot$$ (53)

These can be used to find the effect on the elements of any number of perturbations which are derived from a conservative potential, such as the gravitational interaction potential given by Eq. (13). That is, setting $\mathcal{H}^{(1)}_{\rm r}=\Psi$, we utilize Eqs. (52)-(53) to find that

      

$$\dot{\alpha'_1} -\frac{\mu(x\cos t-2y\sin t)}{r^3}$$ = (54)

$$\dot{\beta'_1} \frac{\mu(x\sin t+2y\cos t)}{r^3}$$ = (55)

$$\dot{\alpha'_2} -\frac{\mu z \cos t}{r^3}$$ = (56)

$$\dot{\beta'_2} \frac{\mu z \sin t}{r^3}$$ = (57)

$$\dot{\alpha'_3} -\frac{\mu y}{r^3}$$ = (58)

$$\dot{\beta'_3} \frac{\mu(2x-3yt)}{r^3}$$ = (59)

where $x=x(\Xi',t),~y=y(\Xi',t),~z=z(\Xi',t),~r=r(\Xi',t)$ are functions of the modified epicyclic elements as defined by Eqs. (49)-(51). Equations (54)-(59) constitute canonical Hill equations. They can be re-written in the state-space form

$$\Sigma:\dot{\Xi'}(t)=F(\Xi',t).$$ (60)

These equations can be studied analytically or numerically in order to detect periodic and quasiperiodic orbits for the coorbital dynamics. The next section presents a numerical procedure used to detect such bounded orbits.

6 Numerical search for bounded relative orbits

We illustrate the use of the newly defined orbital elements by implementing a numerical search procedure aimed at detecting bounded planar solutions for the dynamical system $\Sigma$. To this end, we express the modified epicyclic elements using Fourier series expansions of the form

    

$$\alpha'_1(t)=a_0+\sum_{k=1}^\infty a_k^s\sin(kt)+a_k^c\cos(kt),$$ (61)

$$\alpha'_3(t)=b_0+\sum_{k=1}^\infty b_k^s\sin(kt)+b_k^c\cos(kt),$$ (62)

$$\beta'_1(t)=c_0+\sum_{k=1}^\infty c_k^s\sin(kt)+c_k^c\cos(kt),$$ (63)

$$\tilde{\beta'_3}(t)\doteq\beta'_3(t)-3\alpha'_3(t)t={\rm d}_0+\sum_{k=1}^\infty {\rm d}_k^s\sin(kt)+{\rm d}_k^c\cos(kt),$$ (64)

and perform a numerical search using a gradient-based parameter optimization routine, comprising the following steps. First, we substitute Eqs. (61)-(64), truncated at some given l, into Eqs. (54)-(59). The unknown parameters are:

 

$$\vec a = \{a_0,a_1^s,\ldots a_l^s,a_1^c,\ldots a_l^c,b_0,b_1^s,\ldots b_l^s,b_1^c,\ldots b_l^c,$$ (65)

$$c_0,c_1^s,\ldots c_l^s,b_1^c,\ldots c_l^c~d_0,d_1^s, \ldots d_l^s,d_1^c,\ldots d_l^c\}.$$

Next, we find the optimal Fourier series coefficients by performing the unconstrained least-squares minimization

$$\vec{a}^*= \arg\min_{\vec{a}}\int_0^{t_f}\Vert\dot\Xi'(\vec{a},t)-F([\Xi'(\vec{a},t)]\Vert^2{\rm d}t.$$ (66)

After the optimization routine converges to a (possibly local) minimum, we calculate the initial conditions from:

    

$$\alpha'_1(0) a_0+\sum_{k=1}^l a_k^c$$ = (67)

$$\alpha'_3(0) b_0+\sum_{k=1}^l b_k^c$$ = (68)

$$\beta'_1(0) c_0+\sum_{k=1}^l c_k^c$$ = (69)

$$\beta'_3(0) d_0+\sum_{k=0}^l d_k^c$$ = (70)

and use them to integrate the differential Eqs. (60). If the minimum found is larger than zero, there will be differences between the time histories of the epicyclic elements as obtained from the simulation, and Eqs. (61)-(64). However, after running a sufficient number of random initial guesses, there is a large ensemble of solutions yielding a minimum which is sufficiently close to zero, implying that the integrated solutions and the static solutions match.

Table 1: Initial conditions of modified epicyclic elements.

We emphasize that the optimization procedure described above is static. That is, the differential Eqs. (54)-(59) are transformed into algebraic equations using the pre-defined, periodic, topology of the (candidate) solutions given by Eqs. (61)-(64).

To illustrate the results obtained using the described methodology, we chose $\mu=4.518284\times 10^{-9}$, representing the coorbital system of the Saturnian satellites Epimetheus and Janus, and randomly selected initial guesses for $\vec a^*$. A number of retrograde quasiperiodic orbits were found. Figure 1 depicts some of these orbits. The top plots in this figure describe the orbits in the x-y plane, i.e. the configuration space, and the bottom plots describe the orbits in terms of the guiding center, defined by $(2\alpha'_3, \tilde{\beta}'_3)$. All the axes are normalized by $\epsilon$. The motion is coorbital if the guiding center is contained within the annulus of the colinear Lagrangian equilibrium points L₁ and L₂. In this case, the dynamics are determined by the 1:1 mean motion commensurability.

	Figure 1: Quasiperiodic satellite orbits in the x-y plane (top plots) and their guiding centers in the $\tilde{\beta}'_3-\alpha'_3$ plane (bottom plots).
Open with DEXTER

Orbits (a) represents motion with small magnitude of the guiding center (meaning small variation of the normalized coordinate $2\alpha_3'/\epsilon$), orbit (b) represents motion with medium magnitude, and orbits (c) and (d) represent motion with large magnitude of the guiding center, which is still contained within the annulus of the Lagrangian equilibrium points. Similar result where obtained by Namouni (1999) using classical orbital elements. The initial conditions used to integrate the equations of motion (60) are given in Table 1. We emphasize that by selecting suitable initial conditions, the center of mass of the coorbital satellites will follow the reference unit circle. The analysis of the relative motion is then carried out relative to the known motion of the center of mass.

7 Conclusions

This paper developed a Hamiltonian framework for the analysis of coorbital motion in terms of canonical relative motion elements we termed "epicyclic'' orbital elements. The epicyclic elements are constants of the motion.

We conclude that the approach presented above is very useful for modelling coorbital motion. It renders an analytic insight to the coorbital motion problem and yields a convenient framework for numerical analysis of the coorbital motion problem.

There are many extensions to this approach. We are currently working on developing a solution for coorbital motion about elliptical orbits. We are also looking at the perturbations due to zonal and tesseral gravitational harmonics.

References

 
• Broucke, R. A. 1999, Cel. Mech. Dyn. Astron., 73, 281 In the text NASA ADS  
• Goldstein, H. 1980, Classical Mechanics (New Jersey: Addison-Wesley) In the text  
• Hénon, M. 1969, A&A, 1, 223 In the text NASA ADS  
• Hénon, M., & Petit, J. M. 1986, Celes. Mech., 38, 67 In the text NASA ADS  
• Hill, G. W. 1878, Am. J. Math., 1, 5 In the text  
• Libre, J., & Ollé, M. 2001, A&A, 378, 1087 In the text NASA ADS  
• Namouni, F. 1999, Icarus, 137, 293 In the text NASA ADS  
• Salo, H., & Yoder, C. F. 1988, A&A, 205, 309 In the text NASA ADS  
• Scheeres, D. J., Hsiao, F. Y., & Vinh, N. X. 2003, J. Guidance, Control Dynamics, 26, 62 In the text  
• Villac, B., & Scheeres, D. J. 2003, J. Guidance, Control, Dynamics, 26, 224 In the text  
• Winter, O. C., & Murray C. D. 1997, A&A, 319, 290 In the text NASA ADS  
• Yoder, C. F., Colombo, G., Synnott, S. P., & Yoder K. A. 1983, Icarus, 53, 431 In the text NASA ADS