Hill Problem Analytical Theory to the Order Four: Application to the Computation of Frozen Orbits around Planetary Satellites

Abstract

Frozen orbits of the Hill problem are determined in the double-averaged problem, where short and long-period terms are removed by means of Lie transforms. Due to the perturbation method we use, the initial conditions of corresponding quasi-periodic solutions in the nonaveraged problem are computed straightforwardly. Moreover, the method provides the explicit equations of the transformation that connects the averaged and nonaveraged models. A fourth-order analytical theory is necessary for the accurate computation of quasi-periodic frozen orbits.

1. Introduction

Besides its original application to the motion of the Moon [1], the Hill problem provides a good approximation to the real dynamics of a variety of systems, encompassing the motion of comets, natural and artificial satellites, distant moons of asteroids, or dynamical astronomy applications [2–4]. Specifically, the Hill model and its variations [5–9] are useful for describing the motion about planetary satellites. In addition, the Hill problem is an invariant model that does not depend on any parameter, thus, giving broad generality to the results, whose application to different systems becomes a simple matter of scaling. Note that Hill's case of orbits close to the smaller primary is a simplification of the restricted three-body problem, which in turn is a simplification of real models.

A classical result shows that low eccentricity orbits around a primary body are unstable for moderate and high inclinations due to third-body perturbations [10]. Almost circular orbits close to the central body remain with low eccentricity in the long-term only when the mutual inclination with the perturbing body is less than the critical inclination of the third-body perturbations 𝐼=39.2^∘ (see [11] and references therein). Because of their low eccentricity, high inclination orbits are precisely the candidate orbits for science missions around natural satellites. Therefore, a good understanding of the unstable dynamics of the Hill problem is required.

The study of the long-term dynamics is usually done in the double-averaged problem. After removing the short- and long-period terms, and truncating higher-order terms, the problem is reduced to one degree of freedom in the eccentricity and the argument of the periapsis. As the double-averaged problem is integrable and the corresponding phase space is a compact manifold, the solutions are closed curves and equilibria. The latter are orbits that, on average, have almost constant eccentricity and fixed argument of the periapsis, and are known as frozen orbits.

To each trajectory of the double-averaged problem it corresponds a torus of quasiperiodic solutions in the nonaveraged problem. The accurate computation of initial conditions on the torus requires the recovery of the short- and long-period effects that were eliminated in the averaging. This is normally done by trial and error, making iterative corrections on the orbital elements, although other procedures can be applied [12].

Our analytical theory is computed with Deprit's perturbation technique [13]. The procedure is systematic and has the advantage of providing the explicit transformation equations that connect the averaged analysis with proper initial conditions of the nonaveraged problem. A second-order truncation of the Hamiltonian shows that there are no degenerate equilibria and, therefore, it is sufficient to give the qualitative description of the reduced system. However, the second-order truncation introduces a symmetry between the direct and retrograde orbits that is not part of the original problem, and a third-order truncation is required to reveal the nonsymmetries of the problem.

While, in general, the third-order theory provides good results in the computation of quasiperiodic, frozen orbits, its solutions are slightly affected by long-period oscillations. This fact may adversely affect the long-term evolution of the frozen orbits and it becomes apparent in the computation of science orbits about planetary satellites, a case in which small perturbations are enough for the unstable dynamics to defrost the argument of the periapsis. Then, the orbit immediately migrates along the unstable manifold with an exponential increase in the eccentricity.

We find that a higher-order truncation is desirable if one wants to use the analytical theory for computing accurate initial conditions of frozen orbits. The computation of the fourth-order truncation removes almost all adverse effects from the quasiperiodic solutions, and shows a high degree of agreement between the averaged and nonaveraged models even in the case of unstable orbits.

Whereas the third-body perturbation is the most important effect in destabilizing science orbits around planetary satellites, the impact of the nonsphericity of the central body may be taken into account. The previous research including both effects has been limited up to third-order theories (see [14] and references therein), but from the conclusions of this paper it may worth to develop a higher-order theory including the inhomogeneities of the satellite's gravitational potential.

2. Double-Averaged Hill Problem to the Fourth-Order

The equations of motion of the Hill problem are derived from the Hamiltonian

ℋ=12(𝐗⋅𝐗)−𝝎⋅(𝐱×𝐗)+𝑊(𝐱),𝑊=𝜔2²𝑟²−3𝑥²−𝜇𝑟,(2.1) where, in the standard coordinate system of Hill's model, 𝐱=(𝑥,𝑦,𝑧) is the position vector, 𝐗=(𝑋,𝑌,𝑍) is the vector of conjugate momenta, 𝑟=||𝐱||, and both the rotation rate of the system 𝜔=||𝝎|| and the gravitational parameter 𝜇 of the primary are set to 1 in appropriate units.

The problem is of three degrees of freedom, yet admitting the Jacobi constant ℋ=−𝐶/2. Despite its nonintegrability, approximate solutions that explain the long-term dynamics can be found by perturbation methods. Close to the central body the Hill problem can be written as the perturbed two-body problem

ℋ=12𝑋²+𝑌²+𝑍²−1𝑟−𝜖(𝑥𝑌−𝑦𝑋)+𝜖2²𝑟²−3𝑥²,(2.2) where the first three terms of Hamiltonian (2.2) correspond to the Keplerian motion in the rotating frame and 𝜖 is a formal parameter introduced to manifest the importance of each effect. Thus, the Coriolis term is a first order effect and the third-body perturbation appears at the second-order .

To apply perturbation theory, we formulate the problem in Delaunay variables (ℓ,𝑔,ℎ,𝐿,𝐺,𝐻), where ℓ is the mean anomaly, 𝑔 is the argument of the periapsis, ℎ the argument of the node in the rotating frame, 𝐿=√𝜇𝑎 is the Delaunay action, 𝐺=𝐿√1−𝑒² is the modulus of the angular momentum vector, 𝐻=𝐺cos𝐼 is its polar component, and 𝑎,𝑒,𝐼, are usual orbital elements: semimajor axis, eccentricity, and inclination.

Our theory is based on the use of Lie transforms as described by Deprit [13, 15]. It has the advantage of connecting the averaged and original problems through explicit transformation equations. After removing the short- and long-period terms we get the transformed Hamiltonian

𝒦=𝐾_0,0+𝜀𝐾_0,1+𝜀²2𝐾^0,2+𝜀³6𝐾^0,3+𝜀⁴24𝐾^0,4,(2.3) where 𝜀=𝐿³,

𝐾_0,0=−1𝐾_0,12𝐿²,𝐾=𝐾_0,02𝜎,_0,2=𝐾_0,0142+3𝑒²2−3𝑠²+15𝑒²𝐾𝑠²cos2𝑔_0,3,=𝐾_0,02732𝜎2𝑠²𝐾+50−17𝑠²𝑒²+15𝑒²𝑠²cos2𝑔_0,4,=𝐾_0,0−35123285𝑠⁴+8784−𝑒⁴7c0o8s𝑠4₂𝑔−12𝑠²3996−2940𝑠²−4582−4035𝑠²𝑒²𝑒²cos2𝑔−9𝑠⁴−144926−941𝑠²+244𝑠⁴𝑒²+910728−15208𝑠²+5007𝑠⁴𝑒⁴,(2.4) and 𝜎=𝐻/𝐿=𝑐𝜂,𝜂=√1−𝑒² is the eccentricity function, and we use the abbreviations 𝑠≡sin𝐼,𝑐≡cos𝐼. Details on the perturbation method and expressions to compute the transformation equations of the averaging are given in the appendix.

The double-averaged Hamiltonian (2.3) depends neither on the mean anomaly nor on the argument of the node. Therefore, the corresponding conjugate momenta, 𝐿 and 𝐻, are integrals of the reduced problem and the Hamiltonian (2.3) represents a one degree of freedom problem in 𝑔 and 𝐺. The equations of motion are computed from Hamilton equations d𝑔/d𝑡=𝜕𝒦/𝜕𝐺,d𝐺/d𝑡=−𝜕𝒦/𝜕𝑔,

d𝑔d𝑡=𝐺65𝑐²−𝜂²−5𝑐²−𝜂²cos2𝑔+27𝜀𝜎4𝐺5𝑐²+3𝜀₂+11𝜂²−5𝑐²−𝜂²cos2𝑔128𝐺2113𝑐²−3285𝑐⁴+3915+9165𝑐⁴𝜂²+1581+7791𝑐²𝜂⁴−+41089052𝑒₂𝑐²−1095𝑐⁴+19+2565𝑐⁴𝜂²−547+1744𝑐²𝜂⁴cos2𝑔d𝐺𝑠²𝑐²−𝜂²cos4𝑔,d𝑡=−34𝑒²𝑠²+𝜀5(8+9𝜀𝜎)sin2𝑔₂322509−1095𝑐²+547+4035𝑐²𝜂²sin2𝑔−1095𝑒²𝑠²sin4𝑔.(2.5) Once 𝑔 and 𝐺 are integrated for given initial conditions, the secular variations of ℓ and ℎ are computed from simple quadratures derived from Hamilton equations dℎ/d𝑡=𝜕𝒦/𝜕𝐻, dℓ/d𝑡=𝜕𝒦/𝜕𝐿,

ℎ=ℎ₀+𝜕𝜕𝐻𝒦(𝑔(𝑡),𝐺(𝑡);𝐻,𝐿)d𝑡,ℓ=ℓ⁰+𝜕𝜕𝐿𝒦(𝑔(𝑡),𝐺(𝑡);𝐻,𝐿)d𝑡.(2.6)

3. Qualitative Dynamics

The flow can be integrated from the differential equations mentioned previously, (2.5). However, since the system defined by (2.5) is integrable, the flow is made of closed curves and equilibria, and it can be represented by contour plots of Hamiltonian (2.3). Thus, for given values of the dynamical parameters 𝐿 and 𝐻—or 𝜀 and 𝜎—we can plot the flow in different maps that are function of 𝑔, 𝐺. Figure 1 shows an example in semiequinoctial elements (𝑒cos𝑔,𝑒sin𝑔), where we note a hyperbolic point corresponding to an unstable circular orbit, and two elliptic points corresponding to two stable elliptic orbits with 𝑒=0.2 and periapsis at 𝑔=±𝜋/2, respectively.

Figure 1 

Flow in the doubly reduced phase space.

Delaunay variables are singular for zero eccentricity orbits, where the argument of the periapsis and the mean anomaly are not defined, and for equatorial orbits, where the argument of the node is not defined. Hence, it is common to study the reduced phase space in the variables introduced by Coffey et al. [16], see also [8]:

𝜒₁=𝜂𝑒𝑠cos𝑔,𝜒₂=𝜂𝑒𝑠sin𝑔,𝜒₃=𝜂²−121+𝜎²,(3.1) that define the surface of a sphere

𝜒²¹+𝜒²²+𝜒²³=141−𝜎²²(3.2) of radius 𝑅=(1/2)(1−𝜎²) (the sphere representation misses the case 𝐺=𝐻=0, irrelevant in astrodynamics.)

Then, after dropping constant terms and scaling, Hamiltonian (2.3) writes

𝒦=−12𝜂²−30𝜒²²𝜂²+94𝜀𝜎25−24𝜂²−𝜎²−15𝜒²²+𝜀₂𝜂²643815+9528𝜎²+9𝜎⁴−6343+1709𝜎²𝜂²−1824𝜂⁴+6293−821𝜂²−1470𝜎²𝜒²²𝜂²−3285𝜒⁴²𝜂⁴.(3.3) The flow on the sphere is obtained from Liouville equations 𝜒̇_𝑖={𝜒_𝑖;𝒦}, 𝑖=1,2,3, where the dot means derivative in the new time scale. Then,

𝜒̇₁=316𝜂𝜒²645−8𝜂²+𝜀₂+5𝜎²+72𝜀𝜎5−2𝜂²+5𝜎²643815−1824𝜂⁴+6293−821𝜂₂+9528𝜎²+9𝜎⁴−6𝜂²343+1709𝜎²−1470𝜎²𝜒²𝜂²−3285𝜒²𝜒̇₂𝜂⁴,(3.4)=−316𝜂𝜒¹608𝜀²−320+𝜂3⁴6+0𝜀𝜂²𝜎1−2𝜀8₂+576𝜀𝜎+𝜀²343+1709𝜎²293−1470𝜎²𝜒𝜂²²−1095𝜀²𝜒𝜂²⁴𝜒̇₃,(3.5)=38𝜂𝜒¹𝜒₂320+360𝜀𝜎−𝜀²293−1470𝜎²−821𝜂²−1095𝜒²𝜂²,(3.6) with the constraint 𝜒₁𝜒̇₁+𝜒₂𝜒̇₂+𝜒₃𝜒̇₃=0, derived from (3.2).

Equations (3.4)–(3.6) show that circular orbits (𝜒₁=𝜒₂=0, 𝜒₃=𝑅, the “north” pole of the sphere) are always equilibria. Equations (3.5) and (3.6) vanish when 𝜒₁=0,𝜒₂≠0, but (3.4) vanishes only when

1095𝜀²𝜎⁴−𝜎²320+360𝜀𝜎+𝜀²802+2565𝜎²𝜂²+192−216𝜀𝜎−𝜀²362−35𝜎²𝜂⁶−61𝜀²𝜂⁸=0.(3.7) Equation (3.7) is a polynomial equation of degree 8 in 𝜂, therefore admitting eight roots. Note that, for the accepted values of 𝜀≪1, (3.7) is of the form 𝐴²¹−𝐴²²𝑥+𝐴²³𝑥³−𝐴²⁴𝑥⁴=0 that admits a maximum of three real roots, according to Descartes' rule of signs.

The real roots of (3.7) verified by the dynamical constraint |𝜎|≤𝜂≤1 are also equilibria. The root 𝜂=1 marks a change in the number of equilibria due to a “bifurcation” (𝜂>1 could be a root but not an equilibrium). Then, the number of equilibria changes when crossing the line

𝜀=4−27𝜎−45𝜎³±√5076+1473𝜎²+4730𝜎⁴−27375𝜎⁶423+767𝜎²+1470𝜎⁴(3.8) obtained setting 𝜂=1 in (3.7) that establishes a relation between the dynamical parameters 𝜀 and 𝜎 corresponding to bifurcations of circular orbits. Figure 2 shows that this line defines two regions in the parameters plane with different number of equilibria in phase space. Circular orbits in the outside region of the curve are stable. When crossing the line given by (3.8) the number of real roots of (3.7) with dynamical sense increases such that a pitchfork bifurcation takes place: circular orbits change to unstable and two stable elliptic orbits appear with periapsis, respectively at 𝑔=±𝜋/2, as in the example of Figure 1.

Figure 2 

Regions in the parameters plane with different numbers of equilibria.

Note that the curve given by (3.8) notably modifies the classical inclination limit cos²𝐼>3/5 for circular orbits' stability. However, we cannot extend the practical application of the analytical theory to any value of 𝜀. It is common to limit the validity of the Hill problem approximation to one third of the Hill radius 𝑟_H=3^−1/3. Then 𝜀<(𝑟_H/3)^3/2=1/9, including most of the planetary satellites of interest. Figure 3 shows the bifurcation lines of circular orbits in the validity region of the parameters plane with the values of 𝜀 corresponding to low altitude orbits around different planetary satellites highlighted.

Figure 3 

Bifurcation lines in the parameter plane.

A powerful test for estimating the quality of the analytical theory is to check the degree of agreement of the bifurcation lines of the analytical theory with those computed numerically in the nonaveraged problem. To do that we compute several families of three-dimensional, almost circular, periodic orbits of the Hill (nonaveraged) problem that bifurcate from the family of planar retrograde orbits at different resonances. For variations of the Jacobi constant the almost circular periodic orbits evolve from retrograde to direct orbits through the 180 degrees of inclination. At certain critical points, almost circular orbits change from stable to unstable in a bifurcation phenomenon in which two new elliptic periodic orbits appear. The computation of a variety of these critical points helps in determining stability regions for almost circular orbits [17].

The tests done show that the fourth-order theory gives good results for 𝜀<0.05. As presented in Figure 4, the bifurcation line of retrograde orbits clearly diverges from the line of corresponding critical periodic orbits for higher values of 𝜀, and it may be worth developing a higher-order theory that encompasses also the case of Enceladus.

(a)

(b)

Figure 4 

Comparison between the bifurcation line of circular, averaged orbits (full line), and the curve of critical periodic orbits (dots).

4. Frozen Orbits Computation

Hill's case of orbits close to the smaller primary is a simplification of the restricted three-body problem, which in turn is a simplification of real models. Therefore, the final goal of our theory is not the generation of ephemerides but to help in mission designing for artificial satellites about planetary satellites, where frozen orbits are of major interest.

For given values of the parameters𝜀, 𝜎, determined by the mission, a number of frozen orbits may exist. A circular frozen orbit, either stable or unstable, exists always and the computation of real roots |𝜎|≤𝜂≤1 of (3.7), if any, will provide the eccentricities of the stable elliptic solutions with frozen periapsis at 𝑔=±𝜋/2. To each equilibrium of the doubly reduced phase space it corresponds a torus of quasiperiodic solutions in the original, nonaveraged model. In what follows we present several examples that justify the effort in computing a fourth-order theory to reach the quasiperiodicity condition in the Hill problem.

4.1. Elliptic Frozen Orbits

We choose 𝜀=0.0470573,𝜎=0.422618. If we first try the classical double-averaged solution, the Hamiltonian (2.3) is simplified to 𝐾_0,0+𝜀𝐾_0,1+(𝜀²/2)𝐾_0,2, and the existence of elliptic frozen orbits reduces to the case 𝜎²<3/5, 𝑔=±𝜋/2. The eccentricity of the elliptic frozen solutions is then computed from 𝜂=(5𝜎²/3)^1/4 —obtained by neglecting terms in 𝜀 in (3.7). Thus, for the given values of 𝜀 and 𝜎, and taking into account that we are free to choose the initial values of the averaged angles ℓ, ℎ, we get the orbital elements of the first row of Table 1. The left column of Figure 5 shows the long-term evolution of the instantaneous orbital elements for this case, that we call “classical averaging,” in which we find long-period oscillations of more than four degrees in inclination, more than fifteen in the argument of periapsis, and a variation of±0.06 in the eccentricity.

Table 1 

Initial orbital elements of an elliptic frozen orbit for 𝜀=0.0470573, 𝜎=0.422618.

(a)

(b)

Figure 5 

Long-term evolution of the orbital elements of the elliptic frozen orbit.

When computing a second-order theory with the Lie-Deprit perturbation method we arrive exactly at the classical Hamiltonian obtained by a simple removal of the short-period terms and the classical bifurcation condition that results in the critical inclination of the third-body perturbations 𝐼=39.2^∘ [10, 11]. However, now we have available the transformation equations to recover the short- and long-period effects, although up to the first order only. After undoing the transformation equations we find the orbital elements of the second row of Table 1, where we see that all the elements remain unchanged except for the eccentricity and inclination. The long-term evolution of these elements is presented in the right column of Figure 5, in which we notice a significant reduction in the amplitude of long-period oscillations: 2.5^∘ in inclination, around 10^∘ in the argument of the periapsis, and ±0.04in eccentricity.

The results of the third- and fourth-order theories are presented in the last two rows of Table 1 and in Figure 6. The higher-order corrections drive slight enlargements in the semimajor axis. While both higher-order theories produce impressive improvements, we note a residual long-period oscillation in the elements computed from the third-order theory (left column of Figure 6). On the contrary, the orbital elements of the frozen orbit computed with the fourth-order theory are almost free from long-period oscillations and mainly show the short-period oscillations typical of quasiperiodic orbits.

(a)

(b)

Figure 6 

Long-term evolution of the orbital elements of the elliptic frozen orbit.

4.2. Circular Frozen Orbits

If we choose the same value for 𝜀 but now 𝜎=0.777146, frozen elliptic orbits do not exist any longer and the circular frozen orbit is stable. Both the third and fourth-order theories provide good results, but, again, the third-order theory provides small long-period oscillations in the eccentricity whereas the fourth-order theory leads to a quasiperiodic orbit (see Figure 7).

(a)

(b)

(c)

(d)

Figure 7 

Long-term evolution of the orbital elements of the circular stable frozen orbit. (a) and (c) third-order theory. (b) and (d) fourth-order theory.

For 𝜀=0.0339919 and 𝜎=0.34202 the circular frozen orbit is unstable. Due to the instability, a long-term propagation of the initial conditions from either the third or the fourth theory shows that the orbit escapes following the unstable manifold with exponential increase in the eccentricity. But, as Figure 8 shows, the orbit remains frozen much more time when using the fourth-order theory. A variety of tests performed on science orbits close to Galilean moons Europa and Callisto showed that the fourth-order theory generally improves by 50% the lifetimes reached when using the third-order theory.

(a)

(b)

Figure 8 

Long-term evolution of the orbital elements of the circular, unstable, frozen orbit.

4.3. Fourier Analysis

Alternatively to the temporal analysis mentioned previously, a frequency analysis using the Fast Fourier Transform (FFT) shows how initial conditions obtained from different orders of the analytical theory can be affected of undesired frequencies that defrost the orbital elements.

Thus, Figure 9 shows the FFT analysis of the instantaneous argument of the periapsis of the elliptic orbit in the example mentioned previously. Dots correspond to initial conditions obtained from the double-averaged phase space after a classical analysis—that is equivalent to the second-order analytical theory—and the line corresponds to initial conditions obtained from the fourth-order analytical theory after undoing the transformation. While most of the frequencies match with similar amplitudes, in the magnification of the right plot we clearly appreciate a very low frequency of ∼0.15 cycles/year with a very high amplitude in the classical theory that is almost canceled out with the fourth-order approach. The semiannual frequency remains in both theories because it is intrinsic to the problem. It is due to the third-body perturbation and it cannot be avoided.

(a)

(b)

Figure 9 

(a) FFT analysis of the instantaneous argument of the periapsis of the elliptic solution. (b) Magnification over the low frequencies region.

Figure 10 shows a similar analysis for the instantaneous eccentricity of the stable circular orbit mentioned previously. Now, dots correspond to the fourth-order theory and the line to the third-order one (both after undoing the transformation equations). While the third-order theory provides good results, reducing the amplitude of the undesired frequency to low values, the fourth-order theory practically cancels out that frequency.

(a)

(b)

Figure 10 

(a) FFT analysis of the instantaneous eccentricity of the elliptic solution. (b) Magnification over the low frequencies region.

An FFT analysis of unstable circular orbits has not much sense because of the time scale in which the orbit destabilizes.

5. Conclusions

Frozen orbits computation is a useful procedure in mission designing for artificial satellites. After locating the frozen orbit of interest in a double-averaged problem, usual procedures for computing initial conditions of frozen orbits resort to trial-and-error interactive corrections, or require involved computations. However, the explicit transformation equations between averaged and nonaveraged models can be obtained with analytical theories based on the Lie-Deprit perturbation method, which makes the frozen orbits computations straightforward.

Accurate computations of the initial conditions of frozen, quasiperiodic orbits can be reached with higher-order analytical theories. This way of proceeding should not be undervalued in the computation of science orbits around planetary satellites, a case in which third-body perturbations induce unstable dynamics.

Higher-order analytical theories are a common tool for computing ephemeris among the celestial mechanics community. They are usually developed with specific purpose, sophisticated algebraic manipulators. However, the impressive performances of modern computers and software allow us to build our analytical theory with commercial, general-purpose manipulators, a fact that may challenge aerospace engineers to use the safe, well-known techniques advocated in this paper.

Appendix

Let 𝒯∶(𝐱,𝐗)→(𝐱^,𝐗^), where 𝐱 are coordinates and 𝐗 their conjugate momenta, be a Lie transform from “new” (primes) to “old” variables. If 𝑊=∑_𝑖(𝜖^𝑖/𝑖!)𝑊_𝑖+1(𝐱,𝐗) is its generating function expanded as a power series in a small parameter 𝜖, a function ℱ=∑_𝑖(𝜖^𝑖/𝑖!)𝐹_𝑖,0(𝐱,𝐗) can be expressed in the new variables as the power series (𝒯∶ℱ)=∑_𝑖(𝜖^𝑖/𝑖!)𝐹_0,𝑖(𝐱^,𝐗^) whose coefficients are computed from the recurrence

𝐹_𝑖,𝑗=𝐹_{𝑖+1,𝑗−1}+_{0≤𝑘≤𝑖}𝑘𝑖𝐹^{𝑘,𝑗−1};𝑊_{𝑖+1−𝑘},(A.1) where {𝐹_{𝑘,𝑗−1};𝑊_{𝑖+1−𝑘}}=∇_𝐱𝐹_{𝑘,𝑗−1}⋅∇_𝐗𝑊_{𝑖+1−𝑘}−∇_𝐗𝐹_{𝑘,𝑗−1}⋅∇_𝐱𝑊_{𝑖+1−𝑘}, is the Poisson bracket. Conversely, the coefficients 𝑊_𝑖+1 of the generating function can be computed step by step from (A.1) once corresponding terms 𝐹_0,𝑖 of the transformed function are chosen as desired. In perturbation theory it is common to chose the 𝐹_0,𝑖 as an averaged expression over some variable, but it is not the unique possibility [18]. Full details can be found in the literature [19, 20].

To average the short-period effects we write Hamiltonian (2.2) in Delaunay variables as

ℋ=𝐻_0,0+𝜖𝐻_1,0+𝜖²2𝐻^2,0+𝜖³6𝐻^3,0+𝜖⁴24𝐻^4,0,(A.2) where 𝐻_0,0=−1/(2𝐿²), 𝐻_1,0=−𝐻, 𝐻_2,0=𝑟²{1−3[cos(𝑓+𝑔)cosℎ−𝑐sin(𝑓+𝑔)sinℎ]²}, and 𝐻_3,0=𝐻_4,0=0. Note that the true anomaly 𝑓 is an implicit function of ℓ.

Since the radius 𝑟 never appears in denominators, it results convenient to express Hamiltonian (A.2) as a function of the elliptic—instead of the true—anomaly 𝑢 by using the ellipse relations 𝑟sin𝑓=𝜂𝑎sin𝑢, 𝑟cos𝑓=𝑎(cos𝑢−𝑒), 𝑟=𝑎(1−𝑒cos𝑢).

After applying the Delaunay normalization [21] up to the fourth-order in the Hamiltonian, we get

ℋ^=𝐻_0,0+𝜖𝐻_0,1+𝜖²2𝐻^0,2+𝜖³6𝐻^0,3+𝜖⁴24𝐻^0,4,(A.3) where, omitting primes,

𝐻_0,0=−12𝐿²𝐻_0,1,𝐻=𝐻_0,0𝜀2𝑐𝜂,_0,2=𝐻_0,0𝜀8²4+6𝑒²2−3𝑠²+15𝑒₂+3𝑠²cos2ℎ𝐻2𝑠²cos2𝑔+(1−𝑐)²cos(2𝑔−2ℎ)+(1+𝑐)²cos(2𝑔+2ℎ)_0,3,=𝐻_0,0458𝜀³𝑒²𝐻𝜂(1+𝑐)²cos(2𝑔+2ℎ)−(1−𝑐)²cos(2𝑔−2ℎ)_0,4,=𝐻_0,0×16−45731𝜀+2⁴282𝑐₂+63𝑐⁴−144227+90𝑐²+59𝑐⁴𝑒²+−1284𝑠₂227+610𝑐²−701𝑐⁴𝑒⁴−24𝑠²558+270𝑐²+109−555𝑐²𝑒²𝑒²cos2𝑔−48(1+21𝑐6)₂+56𝑐²−8161+59𝑐²𝑒²−11−701𝑐²𝑒⁴cos2ℎ−48(1−𝑐)₂338−90𝑐+90𝑐²−91−185𝑐+185𝑐²𝑒²𝑒²cos(2𝑔+2ℎ)+6𝑠₄338+90𝑐+90𝑐²−91+185𝑐+185𝑐²𝑒²𝑒²cos(2𝑔−2ℎ)−60𝑠₂56−472𝑒²+701𝑒⁴cos4ℎ+1710𝑠⁴𝑒⁴cos4𝑔+1140𝑠₂18−37𝑒²𝑒²(1+𝑐)²cos(2𝑔+4ℎ)+(1−𝑐)²cos(2𝑔−4ℎ)+285𝑒₄𝑒⁴(1+𝑐)²cos(4𝑔+2ℎ)+(1−𝑐)²cos(4𝑔−2ℎ)(1+𝑐)⁴cos(4𝑔+4ℎ)+(1−𝑐)⁴cos(4𝑔−4ℎ).(A.4) The generating function of the transformation is 𝑊=𝑊₂+(1/2)𝑊₃, where

𝑊₂=𝐿𝜀²×41922−3𝑠₂+6𝑠₂3𝑒5+3𝜂²𝑆_1,0,0−9𝑒²𝑆_2,0,0+𝑒³𝑆_3,0,0+6𝑠₂𝑒35+3𝜂²𝑆_1,0,2+𝑆_1,0−2−9𝑒𝑆_2,0,2+𝑆_2,0−2+𝑒²𝑆_3,0,2+𝑆_3,0−2+6𝑠₂(1+𝜂)²15𝑒𝑆_1,2,0−(9−6𝜂)𝑆_2,2,0+𝑒𝑆_3,2,0+3(1(+1−𝑐)₂𝜂)²15𝑒𝑆_1,−2,0−(9+6𝜂)𝑆_2,−2,0+𝑒𝑆_3,−2,0+3(1−𝑐)₂(1+𝜂)²15𝑒𝑆_1,2,2−(9−6𝜂)𝑆_2,2,2+𝑒𝑆_3,2,2+3(1−𝑐)₂(1+𝜂)²15𝑒𝑆_1,2,−2−(9−6𝜂)𝑆_2,2,−2+𝑒𝑆_3,2,−2+3(1+𝑐)₂(1−𝜂)²15𝑒𝑆_1,−2,2−(9+6𝜂)𝑆_2,−2,2+𝑒𝑆_3,−2,2𝑊(1−𝜂)²15𝑒𝑆_1,−2,−2−(9+6𝜂)𝑆_2,−2,−2+𝑒𝑆_3,−2,−2₃,=𝐿𝜀³×722𝑒5𝑠6₂+88𝑒₃13+3𝜂²𝑆_1,0,2−𝑆_1,0,−2−24𝑒²𝑠²17+4𝜂²𝑆_2,0,2−𝑆_2,0,−2+36𝑒(𝑠1²+𝑆𝜂)_3,0,21−3+𝑆𝜂+_3,0,−28𝜂−₂6𝑒⁴𝑠²𝑆_4,0,2−𝑆_4,0,−2+36𝑒(1−𝜂)13−𝜂+8𝜂₂(1+𝑐)²𝑆_1,2,2−(1−𝑐)²𝑆_1,2,−2−12(1+𝜂)₂(1−𝑐)²𝑆_1,−2,2−(1+𝑐)²𝑆_1,−2,−2−12(1−𝜂)₂17−6𝜂−8𝜂²(1+𝑐)²𝑆_2,2,2−(1−𝑐)²𝑆_2,2,−2+4(1+𝜂)₂17+6𝜂−8𝜂²(1−𝑐)²𝑆_2,−2,2−(1+𝑐)²𝑆_2,−2,−2+4(1−𝜂)₂𝑒(11−6𝜂)(1+𝑐)²𝑆_3,2,2−(1−𝑐)²𝑆_3,2,−2−3(1+𝜂)₂𝑒(11+6𝜂)(1−𝑐)²𝑆_3,−2,2−(1+𝑐)²𝑆_3,−2,−2−3(1−𝜂)₂𝑒²(1+𝑐)²𝑆_4,2,2−(1−𝑐)²𝑆_4,2,−2𝑒²(1−𝑐)²𝑆_4,−2,2−(1+𝑐)²𝑆_4,−2,−2.(A.5) We shorten notation calling 𝑆_{𝑖,𝑗,𝑘}≡sin(𝑖𝑢+𝑗𝑔+𝑘ℎ).

The Lie transform of generating function 𝑊 can be applied to any function of Delaunay variables ℱ=∑_𝑖(𝜖^𝑖/𝑖!)𝐹_𝑖(ℓ^,𝑔^,ℎ^,𝐿^,𝐺^,𝐻^). Since 𝑊₁=0, up to the third-order in the small parameter recurrence (A.1) gives

ℱ=𝐹₀+𝜖²2𝐹⁰;𝑊₂+𝜖³6𝐹⁰;𝑊₃.(A.6) Specifically, this applies to the transformation equations of the Delaunay variables themselves, where 𝐹₀∈(ℓ^,𝑔^,ℎ^,𝐿^,𝐺^,𝐻^) and 𝐹_𝑖≡0 for 𝑖>0.

A new application of the recurrence (A.1) to the Hamiltonian 𝒦=∑_0≤𝑖≤4(𝜖^𝑖/𝑖!)𝐾_𝑖,0, where 𝐾_𝑖,0≡𝐻_0,𝑖 of (A.3), allows to eliminate the node up to the fourth-order, obtaining the double-averaged Hamiltonian (2.3). Note that 𝐾_0,4 corrects previous results in [22]. The generating function of the transformation is 𝑉=𝑉₁+𝜀𝑉₂+(𝜀²/2)𝑉₃, where, omitting double primes,

𝑉₁=𝐿3644+6𝑒²𝑠²𝑉sin2ℎ+5(1+𝑐)²𝑒²sin(2𝑔+2ℎ)−5(1−𝑐)²𝑒²sin(2𝑔−2ℎ)₂,=𝐿−3128𝜂6𝑐2−17𝑒²𝑠²+5(1−𝑐)₂sin2ℎ+5(2−9𝑐)(1+𝑐)²𝑒²sin(2𝑔+2ℎ)𝑉(2+9𝑐)𝑒²sin(2𝑔−2ℎ)₃,=𝐿−9×163𝑠2₂768+2𝑠₄456−104𝑐²−8193+754𝑐²𝑒²+47+7831𝑐²𝑒⁴sin2ℎ−32(12+3𝑐2)+₂416𝑒²−1803𝑒⁴sin4ℎ+32(1−𝑐)₂𝑒²2323−285𝑐+780𝑐²−527−1135𝑐+2125𝑐²𝑒²sin(2𝑔+2ℎ)+220𝑠₂𝑒²2323+285𝑐+780𝑐²−527+1135𝑐+2125𝑐²𝑒²sin(2𝑔−2ℎ)+4520𝑠𝑒²₂4−11𝑒²(1+𝑐)²sin(2𝑔+4ℎ)−(1−𝑐)²sin(2𝑔−4ℎ)−385𝑒₄𝑒⁴(1+𝑐)²sin(4𝑔+2ℎ)−(1−𝑐)²sin(4𝑔−2ℎ)(1+𝑐)⁴sin(4𝑔+4ℎ)−(1−𝑐)⁴sin(4𝑔−4ℎ).(A.7)

The new Lie transform of generating function 𝑉 can be applied to any function of Delaunay variables, and, specifically, to the Delaunay variables themselves. For any 𝜉^∈(ℓ^,𝑔^,ℎ^,𝐿^,𝐺^,𝐻^) the transformation equations of the Lie transform are computed, up to the third-order, from

𝜉^=𝜉^+𝜀𝛿₁+𝜀²2𝛿²+𝜀³6𝛿³,(A.8) where

𝛿₂𝛿₁=𝜉^;𝑉₁,𝛿₃=𝜉^;𝑉₂+𝛿₁;𝑉₁,=𝜉^;𝑉₃+𝜉^;𝑉₂;𝑉₁+𝛿₁;𝑉₂+𝛿₂;𝑉₁.(A.9)

Acknowledgments

This work was supported from Projects ESP 2007-64068 (the first author) and MTM 2008-03818 (the second author) of the Ministry of Science and Innovation of Spain is acknowledged. Part of this work has been presented at 20th International Symposium on Space Flight Dynamics, Annapolis, Maryland, USA, September, 24–28 2007.