A Hopf variables view on the libration points dynamics

TL;DR

This paper analytically investigates the dynamics near libration points in the Hill problem using Lissajous variables and perturbation theory, revealing the structure of invariant manifolds and limitations near libration energies.

Contribution

It introduces a novel analytical approach employing Lissajous variables to simplify and analyze the Hill problem's libration point dynamics.

Findings

Invariant manifolds structure is characterized.

Method reduces problem to a one degree of freedom Hamiltonian.

Computations are accurate near libration point energies.

Abstract

The dynamics about the libration points of the Hill problem is investigated analytically. In particular, the use of Lissajous variables and perturbation theory allows to reduce the problem to a one degree of freedom Hamiltonian depending on two physical parameters. The invariant manifolds structure of the Hill problem is then disclosed, yet accurate computations are limited to energy values close to that of the libration points.

Tables

Table 1. Table 1: Coefficients Q n , j , k subscript 𝑄 𝑛 𝑗 𝑘 Q_{n,j,k} in Eqs. ( 3.3 ) and ( 43 ) .

$Q_{1,0,-2}=ds$	$Q_{1,2,-2}=d^2$                         $\frac{M^N}{1}$
$Q_{1,0,0}=-d^2-s^2$	$Q_{1,2,0}=-2ds$
$Q_{1,0,2}=ds$	$Q_{1,2,2}=s^2$
$Q_{1,1,-3}=\left(7-10{\omega }^2\right)d^{2}s$	$Q_{1,3,-3}=\left(11-2{\omega }^2\right)d^3$
$Q_{1,1,-1}=\left(86-20{\omega }^2\right)d{s}^2+\left(6{\omega }^2+3\right)d^3$	$Q_{1,3,-1}=\left(10{\omega }^2-43\right)d^{2}s$
$Q_{1,1,1}=\left(86-20{\omega }^2\right)d^{2}s+\left(6{\omega }^2+3\right)s^3$	$Q_{1,3,1}=\left(10{\omega }^2-43\right)d{s}^2$
$Q_{1,1,3}=\left(7-10{\omega }^2\right)d{s}^2$	$Q_{1,3,3}=\left(11-2{\omega }^2\right)s^3$

Table 2. Table 2: Coefficients ξ i , j subscript 𝜉 𝑖 𝑗 \xi_{i,j} in Eqs. ( 3.4 ); c 0 = c 3 , 1 − 4 ​ c 3 , 3 subscript 𝑐 0 subscript 𝑐 3 1 4 subscript 𝑐 3 3 c_{0}=c_{3,1}-4c_{3,3} and c 2 = 5 ​ c 3 , 1 − 36 ​ c 3 , 3 subscript 𝑐 2 5 subscript 𝑐 3 1 36 subscript 𝑐 3 3 c_{2}=5c_{3,1}-36c_{3,3} .

$i,j$	$L_{i,j}/\left(4sd\omega \right)$	${\ell }_{i,j}$	$g_{i,j}$	$c_{i,j}$
$1,-3$	$-3c_{1,3}{d}^{2}s$	$c_{1,3}\left(d^2+2s^2\right)d$	$c_{1,3}\left(d^2-2s^2\right)d$ $\frac{M^N}{1}$
$1,-1$	$3c_0{d}^3-2c_{3,1}d{s}^2$	$2c_{3,1}{s}^3-c_{2}s{d}^2$	$8c_{3,1}{d}^{2}s-{\ell }_{1,-1}$
$1,1$	$3c_0{s}^3-2c_{3,1}{d}^{2}s$	$c_{2}d{s}^2-2c_{3,1}{d}^3$	$8c_{3,1}d{s}^2+{\ell }_{1,1}$
$1,3$	$-3c_{1,3}d{s}^2$	$-c_{1,3}\left(2d^2+s^2\right)s$	$c_{1,3}\left(s^2-2d^2\right)s$	$\frac{7}{3}{\omega }^2-\frac{256}{3}$
$2,-1$	$-d^2$	$ds$	$-ds$
$2,0$	$-2ds$	$d^2+s^2$	$d^2-s^2$
$2,1$	$-s^2$	$ds$	$ds$
$3,-3$	$-c_{3,3}{d}^3$	$c_{3,3}{d}^{2}s$	$-c_{3,3}{d}^{2}s$
$3,-1$	$c_{3,1}{d}^{2}s$	$-c_{3,1}\left(d^2+2s^2\right)d$	$c_{3,1}\left(2s^2-d^2\right)d$
$3,1$	$c_{3,1}d{s}^2$	$c_{3,1}\left(2d^2+s^2\right)s$	$c_{3,1}\left(2d^2-s^2\right)s$	$43{\omega }^2-184$
$3,3$	$-c_{3,3}{s}^3$	$-c_{3,3}d{s}^2$	$-c_{3,3}d{s}^2$	$11{\omega }^2-32$

Equations

$$\mathcal{J}=\mbox{$\frac{1}{2}$}(\mbox{\boldmath$P$}\cdot\mbox{\boldmath$P$})-\mbox{\boldmath$N$}\cdot(\mbox{\boldmath$p$}\times\mbox{\boldmath$P$})-\Omega(\mbox{\boldmath$p$}),$$

$$\Omega =\frac{\mu }{R}-\frac{N^2}{2}(R^2-3p_x^2),$$

$$\mathcal{J}=\frac{1}{2}\left(P_x^2+P_y^2+P_z^2\right)+P_x{p}_y-p_x{P}_y-\frac{1}{R}+\frac{1}{2}\left(R^2-3p_x^2\right),$$

$${\dot{p}}_x$$

$${\dot{p}}_y$$

$${\dot{p}}_z$$

$${\dot{P}}_x$$

$${\dot{P}}_y$$

$${\dot{P}}_z$$

$$\rho =3^{-1/3},$$

$${\mathcal{T}}_1:(p_x,p_y,p_z,P_x,P_y,P_z)⟶(x,y,z,X,Y,Z),$$

$$p_x=\rho +x,p_y=y,p_z=z,P_x=X{P}_y=Y+\rho ,P_z=Z.$$

$${\mathcal{T}}_1\circ \mathcal{J}=\mathcal{H}$$

$$\mathcal{H}=\sum _{n\ge 0}\frac{{ϵ}^n}{n!}{H}_n,$$

$$H_{0}=\mbox{$\frac{1}{2}$}(X^{2}+Y^{2})-(xY-Xy)+2(y^{2}-2x^{2})+\mbox{$\frac{1}{2}$}(Z^{2}+4z^{2}),$$

$$H_n=-\frac{n!}{\rho }\frac{r^{n+2}}{{\rho }^{n+2}}{P}_{n+2}(x/r),(n\ge 0),$$

$$\nu =2,$$

$$\left(\begin{array}[]{c}\dot{x}\\ \dot{y}\\ \dot{X}\\ \dot{Y}\end{array}\right)=M_{1}\left(\begin{array}[]{c}{x}\\ {y}\\ {X}\\ {Y}\end{array}\right),\qquad M_{1}=\left(\begin{array}[]{cccccc}0&1&1&0\\ -1&0&0&1\\ 8&0&0&1\\ 0&-4&-1&0\end{array}\right),$$

$${\xi }_i=\sum _{j=1,4}{A}_{i,j}\exp ({\lambda }_{j}t),i=1,\dots 4,$$

$$\lambda =(2\sqrt{7}+1{)}^{1/2},\omega =(2\sqrt{7}-1{)}^{1/2}.$$

$$K_0=\lambda x_1{X}_1+\frac{1}{2}(Y_1^2+{\omega }^2{y}_1^2)+\frac{1}{2}(Z_1^2+{\nu }^2{z}_1^2),$$

$${\mathcal{T}}_2:(x,y,z,X,Y,Z)⟶(x_1,y_1,z_1,X_1,Y_1,Z_1),$$

$$\left(\begin{array}[]{c}x\\ y\\ X\\ Y\end{array}\right)=A\left(\begin{array}[]{c}x_{1}\\ y_{1}\\ X_{1}\\ Y_{1}\end{array}\right),$$

$$A=\left(\begin{array}[]{cccc}2\lambda/\sigma&0&-2\lambda/\sigma&2/\tau\\ (\lambda^{2}-9)/\sigma&-(\omega^{2}+9)/\tau&(\lambda^{2}-9)/\sigma&0\\ (\lambda^{2}+9)/\sigma&(9-\omega^{2})/\tau&(\lambda^{2}+9)/\sigma&0\\ \lambda(\lambda^{2}-7)/\sigma&0&\lambda(7-\lambda^{2})/\sigma&-(\omega^{2}+7)/\tau\end{array}\right),$$

$$J_x=x_1{X}_1,J_y=Y_1^2+{\omega }^2{y}_1^2,J_z=Z_1^2+{\nu }^2{z}_1^2,$$

$${\mathcal{T}}_2\circ \mathcal{H}=\mathcal{K}\equiv \sum _{n\ge 0}\frac{{ϵ}^n}{n!}{K}_n(x_1,y_1,z_1,X_1,Y_1,Z_1),$$

$${\mathcal{T}}_3:(x_1,y_1,z_1,X_1,Y_1,Z_1)⟶(x_2,y_2,z_2,X_2,Y_2,Z_2),$$

$$M_k=Q_k{x}_2^{m_1}{X}_2^{m_2}{y}_2^{m_3}{Y}_2^{m_4}{z}_2^{m_5}{Z}_2^{m_6},k=(m_1,m_2,m_3,m_4,m_5,m_6),$$

$${\mathcal{T}}_3\circ \mathcal{K}=\sum _{n\ge 0}\frac{{ϵ}^n}{n!}{K}_n(y_2,z_2,Y_2,Z_2;J),$$

$$\mathcal{C}=\sum _{n=0}^2\frac{{ϵ}^n}{n!}{C}_n,$$

$$C_0$$

$$C_1$$

$$C_2$$

$$x_1=x_2-{\rho }^2\sigma (\lambda {\Delta }_1+{\Delta }_2),X_1=X_2+{\rho }^2\sigma (\lambda {\Delta }_1-{\Delta }_2),$$

Full text

††journal: Celestial Mechanics and Dynamical Astronomy

1 Introduction

The Hill problem is an useful approximation of the restricted three-body problem —not to be confused with a particular case (Hénon and Petit, 1986)— that, further than its original application to the computation of the orbit of the moon around the earth (Hill, 1878), can be representative of the dynamics of an object under the gravitational atraction of different pairs of solar system bodies. Indeed, when using suitable units of length and time the Hill problem does not depend on any parameter, and, therefore, its application to different scenarios becomes a simple matter of scaling (Szebehely, 1967). In particular, the Hill problem is well suited by itself to study the dynamics about asteroids, yet it may need to be amended to include the important effect of the solar radiation pressure (see García Yárnoz et al., 2015, and references therein). But it can be used too in the description of the most relevant features of the dynamics around planetary satellites (Lara et al., 2007, 2010), a case that may require to further superimpose to the third-body dynamics the nonspherical disturbances of the central body, which can notably modify the orbital behavior close to the origin (see Lidov and Yarskaya, 1974; Vashkovyak, 1996; Scheeres et al., 2001; Lara and San-Juan, 2005; Russell and Lara, 2009, for instance). The Hill problem equations are useful also in the investigation of satellite encounters (Petit and Hénon, 1986), and can capture the bulk of the dynamics of coorbital motion, with different applications to relative spacecraft motion (see Kasdin et al., 2005, and references therein).

The Hill problem has been thoroughly studied numerically by the propagation of periodic and quasi-periodic orbits (Hénon, 1969, 1970, 1974, 2003; Michalodimitrakis, 1980), as well as escape trajectories (Villac and Scheeres, 2003). Other studies provide detailed accounts of the dynamics, including the global description of the planar case for values of the energy corresponding to bounded motion (Simó and Stuchi, 2000). Besides, due to its interest in spacecraft mission design, special emphasis has been given to the study of the stable and unstable manifolds associated to Lissajous orbits, which can be effectively computed from the investigation of the center manifold of the collinear libration points (Gómez et al., 2005; Masdemont, 2005).

The global dynamics of the Hill problem must be necessarily investigated numerically, although purely analytical approaches may provide useful information in those regions of phase space in which the motion can be considered a perturbation of the Keplerian motion —the normalized solution been commonly constrained to the close vicinity of the primary (San-Juan et al., 2006; Lara, 2008; Lara et al., 2010). On the other hand, the normal form approach is not restricted to the case of perturbed Keplerian motion and is customarily used in the computation of the center manifold of the libration points (Gómez et al., 1991). This two degrees of freedom manifold is investigated numerically with the usual tools of non-linear dynamics, as Poincaré surfaces of section and the continuation of periodic and quasi periodic orbits (Gómez et al., 2005). Optionally, analytical approximations to the existing periodic orbits can be obtained with the Lindstedt-Poincaré method (Zagouras and Markellos, 1985) (see, also Farquhar and Kamel, 1973; Richardson, 1980).

Alternatively to the use of Poincaré surfaces of section, the dynamics of the center manifold can be approached analytically, at least for energy values close enough to the energy of the libration points. Indeed, the center manifold Hamiltonian of the Hill problem has the form of a two degrees of freedom perturbed harmonic oscillator in the quasi-resonance condition, which is easily cast into the form of a perturbed elliptic oscillator by the standard introduction of a detuning parameter (Henrard, 1970). Then, the Hamiltonian is rearranged in the form of an unperturbed term in the 1-1 resonance condition, whereas the terms that have been decoupled with the detuning parameter are incorporated into the perturbation. The dynamics of these classes of resonant systems can be efficiently approached analytically using perturbation theory (see the recent review in Marchesiello and Pucacco, 2016) and has been lately applied to the computation of analytical approximations of the Lissajous and Halo orbits in the restricted three-body problem (Celletti et al., 2015). The former was approached by a standard double normalization of the center manifold Hamiltonian in harmonic-type variables, whereas the computation of the later required a preprocessing of the center manifold Hamiltonian in order to apply resonant perturbation theory (Ferraz-Mello, 2007).

On the other hand, the Lissajous transformation (Deprit, 1991) comes out as a convenient option to the customary use of harmonic variables in dealing with elliptic oscillators. It was specifically devised to deal with perturbed elliptic oscillators, and reveals particularly well suited to the construction of a resonant normal form Hamiltonian by standard averaging over the elliptic anomaly. The normalized Hamiltonian is of one degree of freedom and, after reformulation in the Hopf (1931) coordinates, provides a complete description of the reduced dynamics on the sphere (Deprit and Elipe, 1991; Miller, 1991).

In the present research it is shown that application of the Lissajous transformation to the center manifold Hamiltonian of the Hill problem allows for the straightforward construction of an analytical solution for the motion about the libration points. The theory only accommodates the cubic and quartic terms of the perturbation expansion of the Hill problem Hamiltonian about the libration points, leading to an extremely simple normalized Hamiltonian. The standard transformation that performs the reduction of the center manifold is provided up to quadratic corrections, and only short-period terms of comparable accuracy are computed in the normalization of the center manifold Hamiltonian in Lissajous variables. Therefore, the practical application of the current theory is constrained to the lower orders of the energy for which these early truncations made sense. Even so, the insights provided by this simple approach go much further than expected and the solution is able to capture the main features of the dynamics about the libration points. Indeed, it not only shows the existence of the planar and vertical Lyapunov orbits, which exist for all values of the energy above the energy of the libration points; but it also shows the main bifurcations of these fundamental orbits, which occur for energy values considerably far away from that of the libration points. Namely, the bifurcation which gives rise to Halo orbits, and the bifurcation and termination of the two-lane bridge of periodic orbits that connects the families of planar and vertical Lyapunov orbits. The construction of a higher order theory, which will notably improve the accuracy of the solution, is just a matter of mechanizing computations and is not discussed here.

The paper is organized as follows. First of all, basic facts of the Hill problem, including information about the main families of periodic orbits related to the libration point dynamics, are recalled in Section 2. Next, the construction of the perturbation solution for the motion in the vicinity of the libration points, which consists of the reduction to the center manifold and the consequent removal of short period effects, is approached is Section 3. It follows the discussion of the reduced phase space in Section 4, where the equilibria of the reduced dynamics are identified with the main existing families of periodic orbits about the libration points. Finally, some validation tests of the proposed solution are presented in Section 5.

2 Hill problem dynamics

In a rotating frame with velocity $N$ in the $z$ axis direction, the $x$ axis defined by the line joining the primaries, the $y$ axis completing a direct frame, and taking one of the primaries as the origin, the Hill problem is defined by the Hamiltonian

[TABLE]

where \mbox{\boldmathp}\equiv(p_{x},p_{y},p_{z}) is position, its conjugate momentum \mbox{\boldmathP}\equiv(P_{x},P_{y},P_{z}) is velocity in the inertial frame, and $\Omega$ is the potential function

[TABLE]

where R=\|\mbox{\boldmathp}\| and $\mu$ is the gravitational parameter.

After scaling units of length by $\mu^{1/3}$ and time by $1/N$, the Hill problem Hamiltonian is rewritten

[TABLE]

showing that the Hill problem does not depend on any parameter. In these non-dimensional units, the Hamiltonian equations stemming from Eq. (2) are

[TABLE]

where over dots denote derivatives in the rotating frame. From Eqs. (3)–(8) it is immediately apparent that planar motions $p_{z}=P_{z}=0$ exist, as well as the two equilibria $\mathcal{L}_{1,2}=\pm(\rho,0,0,0,\rho,0)$, where

[TABLE]

is called the Hill sphere radius, or Hill radius in short. The two equilibria $\mathcal{L}_{1,2}$ are customarily named libration points. Due to the symmetries of the Hill problem with respect to the plane $x=0$, it is enough to discuss the dynamics about just one of the libration points, say $\mathcal{L}_{1}$.

The study of the linearized dynamics about the libration points shows that even though they are unstable equilibria, periodic motion originates from them in the form of small vertical oscillations through them, and planar oscillations around them (Szebehely, 1967). The so called vertical and planar Lyapunov orbits are then grouped into natural families of periodic orbits which are parameterized by the energy, and are customarily computed by numerical continuation techniques (see Doedel et al., 2003, for instance).

The stability of each periodic orbit is characterized by two parameters, say $s_{1}$ and $s_{2}$, where orbit stability requires that both indices are real numbers with absolute value less than 2, and bifurcations of new families of periodic orbits may happen when any of the indices crosses this level. The graphic representation of the evolution of these indices along the family provides useful information. Thus, the stability curves of the family of planar Lyapunov orbits is depicted in Fig. 1, where the usual scaling $2\,\mathrm{arcsinh}\,{s}_{i}/\mathrm{arcsinh}\,2$ is used rather than $s_{i}$. Because the orbits are planar, one of the indices is related to in-plane perturbations whereas the other is related to out-of-pane perturbations. As shown in the figure, planar Lyapunov orbits are highly unstable, and three crossings of the critical value $|2|$ are observed. The first one occurs close to the energy value \mathcal{H}(\mbox{\boldmathp},\mbox{\boldmathP})\approx-2, where the Halo orbits emerge; the second crossing occurs for \mathcal{H}(\mbox{\boldmathp},\mbox{\boldmathP})\approx-0.6, the beginning of a two-lane bridge of periodic orbits which connect the planar and vertical families of Lyapunov orbits. Finally, when \mathcal{H}(\mbox{\boldmathp},\mbox{\boldmathP})\approx 0 a new family of periodic orbits bifurcates with duplication of the period. More details about periodic orbits and other invariant objects of the Hill problem can be consulted in (Gómez et al., 2005).

3 Perturbation solution

First of all, the origin is translated to the libration point by means of the canonical transformation

[TABLE]

given by

[TABLE]

Then, Eq. (1) is rewritten

[TABLE]

For small values of the distance in Hill units $r=\sqrt{x^{2}+y^{2}+z^{2}}$ compared to the Hill radius $\rho\approx 0.7$, the last summand in Eq. (3) can be replaced by the usual expansion in Legendre polynomials. Then, Eq. (3) is written in the form of a perturbation Hamiltonian

[TABLE]

in which the zeroth order term

[TABLE]

is integrable, $\epsilon$ is a formal small parameter indicating the strength of each term $n$ of the perturbation, and

[TABLE]

where $P_{n}$ is the Legendre polynomial of degree $n$.

3.1 Linearized dynamics

The linear dynamics about the libration points is obtained by truncating the perturbation Hamiltonian (13) to retain only the term $H_{0}$. The last summand of Eq. (14) matches the Hamiltonian of the simple harmonic oscillator with frequency

[TABLE]

and, therefore, shows that, in the linear approximation, the motion in ${z}$ and ${Z}$ decouples from the rest of the flow about the libration point and comprises small oscillations in the $z$ axis direction which, therefore, is a stable direction in all cases. More precisely, the projection of the tangent flow in the $({z},{Z})$ plane is made of ellipses, so that, relative to this plane, the equilibria of the Hill problem are of the center type.

On the other hand, from the Hamilton equations of $H_{0}$,

[TABLE]

is a linear differential system with constant coefficients, whose general solution is made of a linear combination of exponentials

[TABLE]

where $\xi_{i}$ stands for $x$, $y$, $X$, and $Y$, respectively, $(A_{i,j})$ is a 4 $\times$ 4 matrix of arbitrary coefficients, which may be expressed as functions of the initial conditions, and the characteristic exponents $\lambda_{j}$ are the eigenvalues of $M_{1}$. Namely $\lambda_{1,2}=\pm\lambda$, $\lambda_{3,4}=\pm\omega\jmath$, with $\jmath=\sqrt{-1}$ and

[TABLE]

Because $\lambda_{1,2}$ are real, the general solution in Eq. (18) has a hyperbolic component. On the other hand, the characteristic exponents $\lambda_{3,4}$ are pure imaginary, thus giving place to an elliptic or center-type component. Therefore, the equilibrium points of the Hill problem are of the saddle $\times$ center $\times$ center type. Notably, the same dynamical behavior results from the Hamiltonian in separate variables

[TABLE]

which is obtained after the canonical transformation

[TABLE]

given by the linear transformation

[TABLE]

in which

[TABLE]

is obtained based on the eigenvector decomposition of $M_{1}$ in Eq. (17), cf. (Jorba and Masdemont, 1999). Note that Eq. (23) is slightly different from the equivalent matrix in (Gómez et al., 2005), where the last column in Eq. (23) appears multiplied by $\omega$. This is just a consequence of the different form of Eq. (20), which is intentionally chosen in preparation for a following Lissajous transformation, from the zero order Hamiltonian in (Gómez et al., 2005).

3.2 Reduction to the center manifold

The flow derived from Eq. (20) admits the three integrals

[TABLE]

which, for given values of each of them, define corresponding invariant manifolds of the linearized motion. When the motion is constrained to the manifold $J_{x}=0$ then $\dot{x}_{1}=\dot{X}_{1}=0$ and the saddle component is removed. Therefore, the manifold $J_{x}=0$ is of the center $\times$ center type, and, for this reason, is called the center manifold.

The existence of the integral $J_{x}$, and, as a consequence, the center manifold, is not limited to the linear dynamics and can be extended to the nonlinear terms of the transformed Hamiltonian

[TABLE]

which is obtained after applying the transformation defined by Eqs. (22)–(23) to all the summands of Eq. (13).

The procedure for extending the integral $J_{x}$ to the nonlinear terms consist in finding a canonical transformation

[TABLE]

that converts Eq. (24) into a normal form such that, for instance (other possibilities may exist, cf. Gómez et al., 2005), in the new variables all the monomials

[TABLE]

with $m_{1}\neq{m}_{2}$ are removed from the Hamiltonian (Giorgilli et al., 1989). This yields a transformed Hamiltonian

[TABLE]

where $J=x_{2}X_{2}$ is an integral. Equation (26) is in the required normal form, and, therefore, has the center manifold $J=0$.

The transformation $\mathcal{T}_{3}$ is computed by canonical perturbation theory (Deprit, 1969). The construction of the center manifold Hamiltonian is simpler and better understood when using complex variables (Kummer, 1976), because they make trivial the solution of the homological equation of the perturbation method. However, this additional change of variables is not necessary for the low order of the theory presented here, and the computations have been made directly in the subindex 2 variables.

After neglecting terms of $\mathcal{O}(\epsilon^{3})$ and higher, the Hamiltonian of the center manifold is obtained by making $J=0$, in Eq. (26), viz.

[TABLE]

with the summands

[TABLE]

Up to the first order, the transformation is given by $y_{1}=y_{2}$, $z_{1}=z_{2}$, $Y_{1}=Y_{2}$, $Z_{1}=Z_{2}$, and

[TABLE]

with

[TABLE]

whose right members must be evaluated using the variables with subindex 2 for the direct corrections in Eq. (31), and using the variables with subindex 1 for the inverse corrections $x_{2}=x_{1}+\rho^{2}\sigma(\lambda\Delta_{1}+\Delta_{2})$, $X_{2}=X_{1}-\rho^{2}\sigma(\lambda\Delta_{1}-\Delta_{2})$.

3.3 Detuning and Lissajous variables

The unperturbed frequency of the oscillations in the $z$ direction can be written as $\nu=\omega\sqrt{1-\delta}$ where, in view of Eqs. (16) and (19),

[TABLE]

is a “detuning” parameter (Henrard, 1970) that amounts to one tenth of the Hill radius $\rho=3^{-1/3}$, and will be taken as a first order perturbation. Then, the Hamiltonian of the center manifold in Eq. (27) is rearranged in the form

[TABLE]

with $\tilde{C}_{1}\equiv{C}_{1}-\frac{1}{2}\omega^{2}\delta z_{2}^{2}$.

Equation (34) can be viewed as the Hamiltonian of a perturbed elliptic oscillator whose principal part comprises two harmonic oscillators in the 1-1 resonance. Moreover, because the perturbation belongs to the real algebra in the Cartesian variables $(y_{2},z_{2},Y_{2},Z_{2})$, the Hamiltonian (34) is advantageously attacked in Lissajous variables (Deprit, 1991).

The Lissajous transformation

[TABLE]

is defined as

[TABLE]

where $s\equiv{s}(L,G;\omega)$ and $d\equiv{d}(L,G;\omega)$ are the state functions

[TABLE]

The variables in Eqs. (36)–(40) have full geometrical meaning: They define an ellipse in the $y_{2}$-$z_{2}$ plane centered at the origin, whose size and shape are defined by the semi-major axis $a$ and semi-minor axis $|b|$ that are derived from the relations $L=\frac{1}{2}\omega(a^{2}+b^{2})$, $G=\omega{a}b$, with the direction of the semi-minor axis with respect to the axis of ordinates defined by the angle $g$; in this ellipse, the elliptic anomaly $\ell$ is measured from the semi-major axis $b$. An analogous ellipse is defined also in the $Y_{2}$-$Z_{2}$ plane, now with semi-major axis $\omega{a}$ and semi-minor one $\omega|b|$ (see p. 209 and ff. of Deprit (1991) for full details).

The Lissajous transformation is applied to Eq. (34), to give

[TABLE]

where

[TABLE]

and the coefficients $Q_{n,j,k}$, which only depend on the momenta $L$ and $G$ through the state functions $s$ and $d$, are given in Table 1.

3.4 Elimination of the eccentric anomaly

The Hamiltonian in Lissajous variables can be reduced to a one degree of freedom Hamiltonian by means of a new canonical transformation

[TABLE]

such that, after truncation to $\mathcal{O}(\epsilon^{2})$,

[TABLE]

where,

[TABLE]

where primes in functions mean the same functions written in the prime variables, and

[TABLE]

are strictly positive irrational numbers ($k_{1}\approx 0.17$, $k_{2}\approx 0.055$, $k_{3}\approx 0.003$, $k_{4}\approx 0.87$, $k_{0}\approx 6.5\times 10^{-7}$) that have been introduced for abbreviating expressions.

Up to the first order, the short-period corrections of the transformation (44) are of the form

[TABLE]

where $\xi\in(\ell^{\prime},g^{\prime},L^{\prime},G^{\prime})$, and $\beta=0$ for the momenta while $\beta=\pi/2$ in the case of the coordinates. The necessary coefficients are given in Table 2, where primes have been dropped for brevity. Note that $G_{i,j}=(m/n)L_{i,j}$ where $m$ is the coefficient of $g^{\prime}$ and $n$ is the coefficient of $\ell^{\prime}$ in the argument of the trigonometric function factored by $G_{i,j}$, and hence the corresponding coefficients are not provided.

4 The reduced phase space

As a result of the averaging the elliptic anomaly beocmes cyclic in Eq. (45), and, therefore, its conjugate momentum $L^{\prime}$ is an integral of the motion. Thus, the problem has been reduced to a one degree of freedom Hamiltonian in the coordinate $g^{\prime}$ and its conjugate momentum $G^{\prime}$.

The $(g,G)$ chart is a cylindrical map that misses the circular orbits representation. Indeed, when the ellipse’s semi-major and semi-minor axes are equal the angle $g$ is undetermined and the Lissajous transformation is singular. This fact does not invalidates the normalization that has been carried out to eliminate $\ell$. Since the perturbation method used is invariant with respect to canonical transformations (Hori, 1966), the resulting theory remains valid after reformulated in nonsingular variables, which, besides, do not need to be canonical, cf. (Deprit and Rom, 1970).

Therefore, circular orbits are not excluded from the theory, and the reduced phase space is conveniently analyzed in terms of the Hopf coordinates. Thus, the new transformation

[TABLE]

defined by

[TABLE]

which projects the reduced phase space onto the sphere centered at the origin

[TABLE]

is applied to the reduced Hamiltonian $\mathcal{B}$, yielding

[TABLE]

where the new abbreviation

[TABLE]

has been introduced.

In the general case, the Hamiltonian flow associated to Eq. (53) is given by the differential system

[TABLE]

The particular case $L^{\prime}=\tilde{L}=\delta^{*}\omega/k_{2}\approx 1.29839$ yields a differential system that is analogous to the Euler equations for the free rigid body motion, and, therefore, can be integrated in terms of Jacobi elliptic functions. The general solution of Eqs. (54)–(56) is not pursued. If found, it would not provide much insight into the nature of the solution due to the unavoidable use of special functions. However, a lot of qualitative and quantitative information can be obtained from the study of particular solutions, like the equilibria, as well as from the graphic representation of the flow.

4.1 Equilibria

When $I_{2}=I_{3}=0$ Eqs. (54)–(56) vanish. In consequence, the points of the sphere

[TABLE]

are always equilibria of the reduced system. On the other hand, when $I_{2}=0$ but $I_{3}\neq 0$, two new equilibria

[TABLE]

may exist, where $I_{1,\mathrm{H}}$ is the solution of $\dot{I}_{2}=0$ with $I_{3}\neq 0$, namely

[TABLE]

which vanishes for the particular value $\tilde{L}$ mentioned before, in which case $I_{3}=\frac{1}{2}\tilde{L}$. In addition, $I_{1,\mathrm{H}}$ needs to fulfill the geometric condition $|I_{1}|\leq\frac{1}{2}L^{\prime}$, from which

[TABLE]

By replacing $L_{0}$ into Eq. (59) it is shown that the new equilibria bifurcate from $E_{-1}$ when $I_{1,\mathrm{H}}=-\frac{1}{2}L_{0}\approx-0.0384303$.

Analogously, when $I_{3}=0$, $I_{2}\neq 0$, two new equilibria

[TABLE]

may exist, where $I_{1,\mathrm{B}}$ is the solution of $\dot{I}_{3}=0$ with $I_{2}\neq 0$, viz.

[TABLE]

which, again, vanishes for $L^{\prime}=\tilde{L}$, a case in which $I_{2}=\frac{1}{2}\tilde{L}$. Note that $(k_{3}+k_{4})I_{1,\mathrm{H}}=2k_{3}I_{1,\mathrm{B}}$. Again, the geometry of the sphere introduces the constraint

[TABLE]

By replacing $L^{\prime}=L_{1}\approx 1.17113$ and $L^{\prime}=L_{2}\approx 1.45668$ into Eq. (62), the sign taken by $I_{1,\mathrm{B}}$ shows that the bifurcation at $L_{1}$ happens from $E_{-1}$, whereas at $L_{2}$ the bifurcation occurs from $E_{+1}$. Thus, the symmetric equilibria $E_{\pm 3}$ migrate from $E_{-1}$ to $E_{+1}$ from increasing values of $L^{\prime}$, or vice-versa when $L^{\prime}$ decreases from $L_{2}$ to $L_{1}$.

The stability of the equilibria can be computed from the usual linearization of the flow in Eqs. (54)–(56). It shows that $E_{1}$ is of the elliptic type for $L^{\prime}<L_{2}$, and then of the hyperbolic type. Besides, $E_{-1}$ is of the hyperbolic type between $L_{0}$ and $L_{1}$, and elliptic otherwise. $E_{\pm 2}$ are stable when they exist ($L^{\prime}>L_{0}$), whereas $E_{\pm 3}$ are unstable from their bifurcation from $E_{-1}$ at $L^{\prime}=L_{1}$ until they collapse with $E_{1}$ at $L^{\prime}=L_{2}$. This behavior will be confirmed in the next section visualizing the flow.

Besides, it is noted that for small enough values of $L^{\prime}$, and, in consequence of $I_{1}$, $I_{2}$, and $I_{3}$, effects of $\mathcal{O}(L^{\prime 2})$ may be neglected. Then the first order truncation of the Hamiltonian (45) can be taken as representative of the motion, namely $\mathcal{B}\approx{B}_{0}+B_{1}$, which, in the Hopf variables is $\mathcal{I}\approx\omega(1-\frac{1}{4}\delta){L}^{\prime}-\frac{1}{2}\delta\omega I_{1}$. Then, $\dot{I}_{1}\approx 0$, $\dot{I}_{2}\approx\frac{1}{2}\delta\omega{I}_{3}$, $\dot{I}_{3}\approx-\frac{1}{2}\delta\omega{I}_{2}$, which show that $I_{1}$ remains approximately constant along each trajectory on the sphere, whereas $I_{2}$ and $I_{3}$ evolve with the slow frequency $\Omega=\frac{1}{2}\delta\omega$ in circumferences parallel to the plane $I_{1}=0$ of radius $\frac{1}{2}(L^{\prime 2}-4I_{1}^{2})^{1/2}$.

A caveat is in order at this point. While the analytical solution has been constrained to the lower orders in the expansion of the Hamiltonian at the libration point, and, in consequence, to compatible orders of the perturbation theories carried out, the bifurcations predicted by the theory occur at high values of $L^{\prime}$. Hence, this information must be taken just as qualitative and keep in mind that the computed bifurcation values may appreciably change when using a higher order theory.

4.2 Visualizing the flow

On the other hand, the changes in the reduced flow discussed previously can be easily visualized in the sphere without need of integrating Eqs. (54)–(56). Indeed, for a given value of the dynamical parameter $L^{\prime}$, a trajectory in the manifold $\mathcal{I}=h$ is defined by the intersection of the two-dimensional surface defined by the Hamiltonian (53) with the surface defined by the sphere (52). Hence, by elimination of $I_{3}$ between these two equations,

[TABLE]

which provides $I_{2}$ as a function of $I_{1}$ and the pair of dynamical parameters $L^{\prime}$ and, $h$. An analogous elimination of $I_{2}$ yields

[TABLE]

which gives $I_{3}\equiv{I}_{3}(I_{1};L^{\prime},h)$. Therefore, each trajectory $\mathcal{I}(I_{1},I_{2},I_{3};L^{\prime})=h$ can be depicted on the sphere from the simple evaluation of the square roots of Eqs. (64) and (65) in such subset of the interval $I_{1}\in[-\frac{1}{2}L^{\prime},\frac{1}{2}L^{\prime}]$ in which the square roots are real.

The sequence of bifurcations of the flow of the reduced problem presented in Fig. 2 has been depicted using this technique. Two different views of each sphere are shown in the figure, the second one obtained by rotating the first one 180 degrees about the axis $I_{3}$. In spite of the various pairs of spheres correspond to different values of the dynamical parameter $L^{\prime}$, and, in consequence, should have different radius, they are represented with a normalized radius 1 to better appreciate the flow. As shown in Fig. 2, for small values of $L^{\prime}$ the two equilibria $E_{\pm 1}$ are stable and the flow circulates about them (first row of Fig. 2). For increasing values of $L^{\prime}$, the flow distorts about $E_{-1}$, until, eventually, $L_{-1}$ changes to instability in a bifurcation event, and two new equilibria appear in the plane $I_{2}=0$ (second row of Fig. 2). The bifurcated equilibria move on the $I_{1}$-$I_{3}$ meridian towards the $\pm{I}_{3}$ axis, while the flow narrows about the plane $I_{3}=0$ plane (third row of Fig. 2). Eventually, the equilibrium $E_{-1}$ comes back to stability in a new bifurcation, and two new unstable equilibria appear in the plane $I_{3}=0$ (fourth row of Fig. 2), which migrate along the equator of the sphere, the $I_{1}$-$I_{2}$ circumference, for increasing values of $L^{\prime}$ until merging in a new bifurcation phenomenon with the equilibrium $E_{+1}$, which undergoes a concomitant change to instability (fifth row of Fig. 2). Further increases of $L^{\prime}$ do not introduce qualitative changes in the flow.

4.3 Orbits of the center manifold

The orbits of the reduced phase space are now identified with orbits of the original Hill problem. From Eqs. (57), (58), and (61), and the inverse transformation of Eq. (51) one finds that $G^{\prime}=0$ for $E_{\pm 1}$ and $E_{\pm 3}$. Hence, except for the short-period effects due to the corrections in Eq. (3.4), the corresponding motion in the center manifold is rectilinear. On the contrary, $G^{\prime}\neq 0$ for $E_{\pm 2}$, which, therefore, correspond to elliptic motion, on average, in the center manifold. Besides, $g^{\prime}=0$ for $E_{1}$, which in view of Eqs. (36)–(39), yields, on average, harmonic oscillations in the $z_{2}$ direction, whereas $g^{\prime}=\frac{\pi}{2}$ for $E_{-1}$ thus constraining the oscillations to the $y_{2}$ direction. On the other hand, the solutions $E_{\pm 2}$ yield elliptic oscillations in the $y_{2}$-$z_{2}$ plane, the area of the respective ellipses depending on $|G^{\prime}|=(L^{\prime 2}-4I_{1,\mathrm{H}}^{2})^{1/2}$. Finally, the equilibria $E_{\pm 3}$ result in rectilinear oscillations in the $y_{2}$-$z_{2}$ plane with inclination $g^{\prime}$ given by the components $I_{1}$ and $I_{2}$ of the equilibria.

When the Cartesian coordinates are recovered through the linear transformation in Eq. (22), one recognizes that the the equilibria of the reduced problem correspond to the well known periodic orbits of the center manifold of the Hill problem (Gómez et al., 2005). Namely,

vertical ($E_{1}$) and planar ($E_{-1}$) Lyapunov orbits

• 2.

Halo orbits, with the upper part towards the libration point ($E_{2}$) and the symmetric one ($E_{-2}$) with the upper part towards the primary

• 3.

orbits of the two lane bridge ($E_{\pm 3}$) connecting vertical and planar Lyapunov orbits

The periods of these periodic orbits is estimated from the rate of variation of the eccentric anomaly in the Lissajous normalized coordinates $T=2\pi/\dot{\ell}^{\prime}$, where $\dot{\ell}^{\prime}=\partial\mathcal{B}/\partial{L}^{\prime}$ is computed from Eq. (45). Up to the first order in the small parameter

[TABLE]

where in each case the values of $g^{\prime}$, $G^{\prime}$ and $L^{\prime}$ must be replaced from those of the corresponding equilibria.

5 Validation tests

Some examples are provided to illustrate the application and performance of the analytical theory.

Thus, starting from $L^{\prime}=0.001$, the Hopf coordinates of the $E_{1}$ equilibrium, corresponding to the Lyapunov vertical solution, are computed from Eq. (57): $I_{1}=0.0005$, $I_{2}=I_{3}=0$, what result into $g^{\prime}=0$ and $G^{\prime}=0$, to which correspond a period $T=3.13965$ of the rectilinear oscillations. Then, for each value $\ell^{\prime}\in[0,2\pi)$, the original Lissajous variables are computed by recovering the short-period corrections in Eq. (3.4). Then, the Lissajous transformation in Eqs. (36)–(39) provides the coordinates $y_{2}$, $z_{2}$, $Y_{2}$, $Z_{2}$, in the center manifold. Figure 3 illustrates how this orbit looks like in the center manifold.

To recover the orbit in the original space, the corrections in Eqs. (3.2) and (33) are computed first to get the subindex 1 variables, and then Eq. (22) provides the Cartesian coordinates of the orbit in the original space, which is illustrated in the left plot of Fig. 4. The analytical solution obtained in this way is obviously periodic by construction (blue points in the left plot of Fig. 4). However, when the initial conditions obtained from the analytical solution for, say, $\ell^{\prime}=0$ are propagated in the original equations of motion, Eqs. (3)–(8), the orbit is not exactly periodic due to the neglected higher order corrections in the transformations computed by perturbation theory (black curve in the left plot of Fig. 4), namely, the reduction to the center manifold and the averaging of the elliptic anomaly, as well as in the computation of the period. Indeed, the differences between the initial state and the computed state after the period given by the analytical approximation is of the order of $10^{-5}$ in this example.

As expected, because the analytical solution is constrained to small values of the distance when compared to the Hill radius, this lack of periodicity is more evident for higher values of $L^{\prime}$, because of the corresponding larger size of the orbit. The differences between the analytical solution and the true orbit for the same initial conditions soon become evident, as illustrated in the right plot of Fig. 4 for $L^{\prime}=0.02$, where dots correspond to the analytical orbit and the black curve represents to the numerical solution. Still, for moderate values of $L^{\prime}$ the analytical solution provides good enough seeds to feed a differential corrections procedure that easily gets the periodic solution. Thus, even though the periodicity error when using initial conditions of the analytical theory is of the order of $10^{-2}$ for $L^{\prime}=0.02$, the differential corrections algorithm in (Lara and Peláez, 2002) only needs to compute four consecutive corrections to converge to a periodic orbit with a periodicity error of $\mathcal{O}(10^{-13})$.

Orbits of the planar Lyapunov family are analytically computed analogously, now starting from the $E_{-1}$ equilibrium in the Hopf coordinates. These orbits get much closer to the primary than corresponding vertical ones for the same values of $L^{\prime}$, and hence the effects of the perturbation are stronger and may manifest the lacks of using a lower order theory much sooner. Therefore, a wise selection of the initial conditions to propagate in the real model may be crucial to obtaining a good approximation of a planar Lyapunov orbit. As shown in Fig. 5, if the initial conditions are taken from the analytical solution for $\ell^{\prime}=0$ (left plot) the highly unstable behavior makes that the orbit very soon departs from the nominal trajectory. On the contrary, choosing $\ell=\pi$ provides a much better approximation of a periodic orbit (right plot), which, again, is easily improved by differential corrections in the original problem.

The situation becomes more critic for increasing values of $L^{\prime}$. Thus, while the shape of the orbit predicted by the analytical solution may resemble the aspect of the corresponding periodic orbit of the Hill problem, the initial conditions obtained from the analytical solution may fail when used as a seed in the search for a true periodic orbit. Because of the lower order of the current theory, this is exactly what happens to Halo orbits derived from the analytical solution, which the actual truncation of the theory to the second order predicts to exist only for $L^{\prime}>0.077$, cf. Eq. (60). As shown in the left plot of Fig. 6, initial conditions of the analytical theory (blue points) clearly fail in providing a true Halo orbit when numerically propagated (black curve), and the situation is even worse for the hodograph (not shown). Still, the size and shape of the analytical orbit is representative of the real Halo trajectory as illustrated in the right plot of Fig. 6, where the black curve is a true periodic orbit of the Halo family that has been computed numerically.

The situation will definitely improve when extending the analytical solution to higher orders. Nevertheless, in view of the periodic orbits of the bridge family linking vertical and planar Lyapunov orbits only exist for much higher values of $L^{\prime}$, cf. Eq. (63), and, therefore, can take values of the $x$ coordinate close to the Hill radius, one should not put big expectations in seizing the real dynamics of these large and highly unstable periodic orbits even when using a very high order theory.

6 Conclusions

A chain of five consecutive canonical transformations (a translation of the origin to the libration point, a decoupling of the linear dynamics, a polynomial normal form, the Lissajous transformation, and a short-period averaging of the elliptic anomaly) reduces the Hamiltonian of he Hill problem to a one degree of freedom Hamiltonian that comprises the dynamics in the vicinity of the libration points. The reduced phase space is advantageously described in the Hopf coordinates, which project the Lissajous variables onto the sphere and provide a deeper insight than the usual surface of section representation. In particular, the change to instability of the Lyapunov planar orbits, with the consequent bifurcation of the family of Halo orbits, is clearly visualized in the Hopf coordinates, and their associated stable and unstable manifolds are unambiguously distinguished. On the contrary, these stable and unstable manifolds almost coincide in the usual surface of section representation, in which the Lyapunov planar orbit is not represented by a point, but by the curve bounding the surface of section, and hence it is difficult to differentiate one manifold from the other.

The reduced dynamics is made of fixed points of the elliptic and hyperbolic types as well as closed curves surrounding these equilibria, which are identified with the well known trajectories of the center manifold. In particular, both types of Lyapunov periodic orbits (planar and vertical) correspond to oscillations that on average remain rectilinear in Lissajous variables, whereas the Lissajous quasi-periodic orbits are, on average, rotating ellipses with variable angular momentum in Lissajous variables.

The accuracy of the analytical solution is limited to energy values close to that of the libration points because of the early truncation of the developments to the fourth order in the ratio distance over Hill radius. However, the range of applicability of the solution could be trivially enlarged by the straightforward computation of higher orders of the perturbation theory.

Acknowledgemnts

The author acknowledges partial support by the Ministry of Economic Affairs and Competitiveness of Spain, under grants ESP2013-41634-P, ESP2014-57071-R and ESP2016-76585-R.