COMPUTING LONG-LIFETIME SCIENCE ORBITS AROUND NATURAL SATELLITES

. Science missions around natural satellites require low eccentricity and high inclination orbits. These orbits are unstable because of the planetary perturbations, making control necessary to reach the required mission lifetime. Dynamical systems theory helps in improving lifetimes reducing fuel consumption. After a double averaging of the 3-DOF model, the initial conditions are chosen so that the orbit follows the stable-unstable manifold path of an equilibria of the 1-DOF reduced problem. Corresponding initial conditions in the non-averaged problem are easily computed from the explicit transformations provided by the Lie-Deprit perturbation method.


1.
Introduction. Mission designing for artificial satellite missions is a two step procedure. First, a deep knowledge of the dynamics is obtained through simplified models, yet retaining the main characteristics of the real problem. The preliminary design in the simplified model is then optimized in a real, ephemeris model.
In the case of planetary satellites, three-body models have been proposed as soon as in the dawn of the space era. Thus, for a Lunar orbiter Kozai [1] deals with a Hill-type problem that includes the non-sphericity of the central body. The recent interest in science missions to some icy moons, as Europa or Enceladus, has resuscitated Kozai's model showing that it accurately fits to the dynamics of other planetary satellites [2,3,4].
Future missions to planetary satellites will require science orbits with low altitude and eccentricity, and high inclination. However, third body perturbations produce unstable behavior of high-inclination orbits around planetary satellites. The instability is manifested by an exponential growth in the orbit eccentricity that soon leads to impact with the surface of the satellite. Therefore, space mission planners must confront the problem of developing strategies to maximize the lifetime of the orbiter -the time after which the orbiter escapes or impacts the planetary satellite.
Dynamical systems theory has proven successful in the past for dealing with unstable motion [5]. The use of unstable periodic orbits and their associated stable and unstable manifolds provides a non-expensive, thrust-free way of designing maneuvers. Also, one can deal with the doubly averaged problem and use the invariant manifolds of selected unstable equilibria to maximize orbit lifetime [6,7].
Classical averaging procedures provide a good understanding of the long term, overall dynamics [2]. However, they lack of providing the equations for transforming the new variables into the old ones and vice-versa, making necessary to perform involved computations to reach this goal (see [6]). On the contrary, the Lie-Deprit perturbation method [8] used in this paper, not only enables the computation of the secular elements to high orders in a systematic way, but it generates explicitly the canonical transformations from the averaged phase space to the original one and its inverses. This is of major relevance to the design of science orbits around planetary satellites because the candidates for nominal orbits obtained from the averaged analysis are easily connected with proper initial conditions for the mission.
The dynamical model of this paper is a Hill-type, invariant problem. From equilibrium theory, it is assumed that the satellite is in synchronous rotation with its orbital motion around the planet and, therefore, the satellite exhibits certain oblateness and eccentricity. The Hamiltonian also considers the effect of a possible latitudinal asymmetry -observed in many solar system bodies-that has relevant effects in the dynamics.
Close to the satellite the ratio rotation rate of the satellite mean motion of the orbiter is a small quantity that naturally scales the Hamiltonian. A double average over the mean anomaly and the argument of the node reduces the problem to one degree of freedom in the eccentricity and the argument of the periapsis, where the equilibria correspond to what aerospace engineers call frozen orbits (see [9] and references therein). The time evolution of the secular terms is obtained to order three in the small parameter. The computation of the initial conditions of long lifetime science orbits from the (second order) transformation equations is exemplified for a Europa orbiter. This paper extends the results of [7] up to the second order. Contrary to [7], we do not start from a Hamiltonian averaged over the short period terms, and the transformation equations of the Delaunay normalization are, therefore, required. Further, the model used in [3] is completed with the effect of the latitudinal asymmetry, but the discussion of the phase space of the new model is deferred to a future paper. Finally, the manifolds of the 3-DOF problem obtained from the analytical theory of this paper should be related to the periodic orbits computed in [4]. The problem of relating the equilibria of the 1-DOF problem with periodic orbits of the non averaged Hamiltonian requires an estimation of the orbit period in the rotating frame and remains open at this stage.

2.
Equations of motion. When considering a synchronously orbiting, planetary satellite whose equator coincides with its orbital plane the problem is invariant in a rotating frame of reference with origin at the center of mass of the satellite. Using Hamiltonian formulation we write where x = (x, y, z) is the position vector of the orbiter, r = ||x||, X = (X, Y, Z) is the vector of conjugated momenta -velocity in the inertial frame, ω is the angular velocity of the satellite around the planet, µ is the satellite's gravitational parameter, and the perturbing function R includes the third body and non-sphericity perturbations. Hamilton equationsẋ = ∇ X H,Ẋ = −∇ x H, of Eq. (1) givë equations of motion corresponding to a nonlinear dynamical system of three degrees of freedom, yet accepting the Jacobi integral H(X, x) = E. We consider the perturbing function where the non-dimensional part N of the perturbing potential is the scaling factor α is the equatorial radius of the satellite, ω = ||ω||, and the parameters J 2 , C 2,2 , and J 3 are the oblateness, ellipticity, and latitudinal asymmetry coefficients respectively, the projective coordinates are r and h is the argument of the node in the rotating frame, I is inclination, θ = f + g, g is the argument of the periapsis, and f is the true anomaly. From equilibrium theory [10], we set C 2,2 = (3/10) J 2 . Finally, we assume J 3 to be of higher order than J 2 .
Hamiltonian (1) is expressed as a perturbed two-body problem where the formal parameter ǫ manifests the order of each perturbation. The zero order is the Keplerian where a is the instantaneous semimajor axis. The first order perturbation corresponds to the Coriolis force where e, is eccentricity, and n is the mean motion of the orbiter. The effects of the oblateness and ellipticity of the satellite, and the planetary perturbation remain at second order of (ω/n) And the effect of the latitudinal asymmetry of the satellite is of third order in (ω/n) For brevity, we use the abbreviations c ≡ cos I, s ≡ sin I.
2 , express the relative importance of the non-sphericity and third-body perturbations. We limit our theory to those regions in phase space where both β and γ are of order 1.
3. Perturbation theory. For our perturbation theory we use Delaunay-type variables (ℓ, g, h, L, G, H) where ℓ is the mean anomaly and the momenta are functions of the orbital elements L = √ µ a, G = L 1 − e 2 , H = G cos I.

MARTIN LARA AND SEBASTIÁN FERRER
We use the Lie transforms perturbation method as described in [8]. Thus, the average over the mean anomaly is performed by a Lie transform which allows to compute the new Hamiltonian in new variables that after truncation to order k is free from the averaged variable (the dash remarks the absence of the corresponding Delaunay element in a function). The Lie-Deprit perturbation methods allows the computation, order by order, of both the averaged Hamiltonian and the required generating function. Thus, to compute the new Hamiltonian averaged in ℓ ′ up to the order 3, we start the procedure by setting Then, we use Deprit's recurrence is the Poisson bracket in Delaunay elements. Full details can be found in [8] (see also [11] where Deprit dealt with normalisations of perturbed Kepler Hamiltonians).
In what follows we omit the "prime" notation without risk of confusion.
After averaging the mean anomaly we obtain where η = G/L is the eccentricity function, and e 2 = 1 − η 2 . Note that β and γ are constants in the 2-DOF problem. At this stage, the polylogarithm function appears in the third order generating function. Until we train ourselves in manipulating this kind of special functions [12] we stop at the third order, meaning that the transformation equations are limited to the second order in the small parameter.
A new Lie transform T 2 : (ℓ ′ , g ′ , h ′ , L ′ , G ′ , H ′ ) → (ℓ ′′ , g ′′ , h ′′ , L ′′ , G ′′ , H ′′ ) eliminates the argument of the node proceeding analogously to the averaging over the mean anomaly, and reduces the 2-DOF system to a 1-DOF. In this way we get the doubly averaged Hamiltonian where the dot means derivative with respect to time. Equations (10)-(11) correspond to an integrable dynamical system of 1-DOF. After finding the solution g(t), G(t), the quadratures providing the secular variation of ℓ and h are obtained from the corresponding Hamilton equationṡ ℓ = ∂K(g(t), G(t)) ∂L ,ḣ = ∂K(g(t), G(t)) ∂H Note that the discussion of the reduced phase space can be done after dropping constant terms in the Hamiltonian (9), and using a new time scale. Then, after scaling, the evolution of the variables (g, η) can be investigated in the four dimensional parametric space (ǫ, β, γ, σ) (see [3] for a detailed discussion for the triparametric case γ ≡ 0).
Delaunay variables are singular for e = 0. Further, the (g, G) representation corresponds to a cylindrical map that presents singularities for circular orbits, and the flow must be properly represented on a closed surface [13]. When rectilinear solutions (η = 0) can be ignored -the case in Astrodynamics, where high eccentricity orbits will impact the satellite-the flow is diffeomorphic to a 2-sphere [14], and it can be analized with Coffey et al. [15] variables χ 1 = η e s cos g, χ 2 = η e s sin g, on the sphere of radius R = χ 2 1 + χ 2 2 + χ 2 3 = 1 2 1 − σ 2 . We delay the discussion of the reduced phase space to a future paper. We just mention that the equilibria of the reduced system Eqs. (10)-(11) correspond to orbits that, on average, have constant eccentricity and argument of periapsis. These orbits are of interest in mission designing for artificial satellites and aerospace engineers usually call them frozen orbits [9]. In general, the frozen orbit values are computed numerically solving the equilibria equationsġ =Ġ = 0 from Eqs. (10)- (11). However, in specific cases frozen orbits can be determined generically.
In this case, we note in Eq. (11) that g = ±π/2 ⇒Ġ = 0. Then, for g = ±π/2 andġ = 0, Eq. (10) establishes a dependence between the eccentricity and the inclination F (e, I; ε, β, γ) = 0 that depends on the parameters [16]. Then, one can follow the evolution of the frozen orbits in the (e, I) plane, as illustrated in Fig. 1, which helps in determining the most convenient solution for different purposes. Note in the example of Fig. 1 that high inclination frozen orbits, say 65 • < I < 105 • , enjoy low eccentricity, say e ∼ 10 −3 , that peaks to almost circular orbits very close to the critical inclination values cos 2 I = 1/5.

5.
Computing long lifetime orbits. An important application of the analytical theory of this paper is the computation of long-lifetime orbits. Low altitude, high inclination orbits around planetary satellites are known to be unstable because of the planetary perturbation. Despite the fact that stability is not discussed in this paper, the graphic representation of the Hamiltonian flow for given values of the parameters can be used to corroborate this fact. Figure 2 provides an example where the parameters were set to β ≈ 0.8, ε ≈ 0.025 (γ = 0) -close to the values of an orbiter 100 km above the surface of the Galilean moon Europa. We note the hyperbolic character of the origin, which corresponds to a circular orbit, with its stable-unstable manifolds highlighted. We also note other equilibria of the flow either stable or unstable: the elliptic points at e cos g = ±0.7, e sin g = 0 (stable frozen orbits with the argument of periapsis g = 0, π, and e = 0.7), and e cos g = 0, e sin g = ±0.89 (stable frozen orbits with the periapsis at g = ±π/2 and e = 0.89), and the hyperbolic points at e cos g = ±0.87, e sin g = 0 (stable frozen orbits with g = 0, π and e = 0.87). We note that the highly elliptic frozen orbits of this example (centers with e = 0.7, e = 0.89, and saddle with e = 0.89) are orbits that impact the satellite. The instability of the frozen orbit of the science mission (corresponding to the origin of Fig. 2) produces that it soon enters the unstable manifold, with a fast increase in the eccentricity that leads to the impact with the surface of the satellite. For a Europa orbiter these instabilities result in relatively short lifetimes, on the order of one month [17]. To enlarge the lifetime of the orbit, one find assistance in Dynamical Systems Theory, and stable-unstable manifolds tours have been proposed [6]. Thus, if initial conditions of a nonimpact orbit are chosen on the stable manifold of the unstable equilibrium, the natural dynamics heads the orbit towards the unstable equilibrium closely over the stable manifold, and then the unstable manifold destabilizes the orbit.
Thus, for given values a, I, we compute the parameters ǫ, β and γ of the desired planet-satellite system. Then, Eqs. (10)-(11) are solved for the equilibria and the lower eccentricity e 0 is selected (note that e 0 = 0 if γ ≡ 0). The averaged Hamiltonian (9) provides the energy k 0 = K(e 0 , g; a, I) of the unstable, frozen orbit, which is also the energy of their associated invariant manifolds. The contour level K(e, g; a, I) = k 0 is, then, the homoclinic trajectory connected with the frozen orbit (dashed lines in Fig. 2).
The maximum eccentricity on a non-impact orbit on the stable manifold of the nominal frozen orbit occurs when the periapsis distance equals the radius of the satellite. This is e M = 1 − (α/a). Once selected e < e M , Eq. (9) K(e M , g; a, I) = k 0 can be solved for g, thus providing the initial conditions of an orbit on the stable manifold of the averaged problem. Note, however, that the initial conditions selected from the averaged problem will not produce, in general, a long lifetime orbit: The short and long period effects missed in the averaging will force the 3-DOF orbit to separate from the averaged manifolds.
The evolution of an orbit whose initial conditions are chosen directly from the averaged dynamics is illustrated in Fig. 3. For the example we use the physical values of Europa, namely µ = 3202.7 km 3 /s 2 , ω = 2.0477 · 10 −5 rad/s, α = 1565 km, and J 2 = 4.355·10 −4 . We choose a = α+100 km, I(e = 0) = 80 degrees; assuming a minimum perigee 2 km larger than the Europa radius, we obtain 1 − (α/a) = 0.06, and set e M = 0.059 and compute the argument of the periapsis from Eq. (9) obtaining g = 33.85 degrees. Finally, we choose h = ℓ = 0 and propagate the initial conditions in the non averaged problem Eq. (2). To guarantee the accurate propagation of the initial conditions we decided to integrate the equations of motion by Recurrent Power Series, and checked that the Jacobi integral E was preserved during the long term propagations to more than 13 significant digits. In the left plot of Fig. 3 we present the radius evolution, where we see that the orbit impacts the surface of Europa after around two months; the plot in the right of Fig. 3 shows that the orbit roughly follows the stable-unstable manifolds. The initial conditions of an orbit of the non-averaged problem, Eq. (1), that, on average, moves on the stable manifold of the frozen orbit, are obtained from the transformation equations provided by the Lie-Deprit perturbation method. Note that we are free to choose h and ℓ, however once the values are selected the corresponding ones in the non averaged problem must be computed. After applying the computed corrections to the orbital elements: ∆a = 1.37 km, ∆e = 0.0011, ∆I = 0.9 • , ∆g = 1 • , the lifetime of the orbit almost doubles the original one. This is appreciated on the left plot of Fig. 4, while the right plot of the Figure shows that, now, on average the orbit moves tightly over the invariant manifolds. Due to the symmetries of the case J 3 = 0, the fact that the path changed from the previous case is irrelevant, because identical lifetimes are obtained when moving over any of the manifolds. Note, however, that it is not true when J 3 = 0, where a bias on the initial argument of the periapsis may be introduced for reaching the longest lifetime [6]. The propagation in 3-DOF of initial conditions of the 1-DOF problem, as in Fig. 3, deserves a final comment. For J 3 = 0, one might find longer lifetimes when g is computed from a second order truncation that when it is computed from the full Hamiltonian (1). We remark that it does not mean that a second order approach provides results better than the third order theory of tis paper. Quite on the contrary, after undoing the Lie transformation equations of the second order approach one finds negligible improvements in the lifetimes. This effect should be attributed to the intricacies of the manifolds of the non averaged model [18] and never to an inaccuracy of our theory. 1 6. Conclusions. Analytical theories are of great use in mission designing for artificial satellites. Modern perturbation methods not only describe the evolution of the secular terms, but provide the explicit transformations that recover the short term dynamics and connect the averaged analysis with proper initial conditions for the mission. This is of major relevance to the design of a science orbit about planetary satellites, because the candidates for nominal orbits live among the set of unstable orbits.
We reach the third order in closed form with our theory. At this stage, special functions appear that need further study to continue the procedure to higher orders if the closed form is desired. However, the third order showed useful in the computation of long lifetime orbits: a simple evaluation of the explicit transformation equations (that are computed literally and once for all) doubles the lifetimes obtained directly from the averaged equations. The Lie-Deprit perturbation method is straightforward and, with the impressive performances of modern computers, it can be easily implemented to higher orders with general purpose commercial algebraic manipulators, as we did in the present work.