Open AccessFull-Length Papers

Local Orbital Elements for the Circular Restricted Three-Body Problem

Luke T. Peterson

https://orcid.org/0000-0003-2648-8094

University of Colorado Boulder, Boulder, Colorado 80309

*Ph.D. Student, Aerospace Engineering Sciences; . Student Member AIAA.

Search for more papers by this author

and

Daniel J. Scheeres

University of Colorado Boulder, Boulder, Colorado 80309

†A. Richard Seebass Endowed Chair Professor, Aerospace Engineering Sciences. Fellow AIAA.

Search for more papers by this author

Published Online:8 Aug 2023https://doi.org/10.2514/1.G007435

Abstract

The Keplerian orbital elements make up a set of parameters that uniquely describe a two-body trajectory. Once perturbations are imposed on two-body dynamics, one often studies the time evolution of the orbital elements. Moving in complexity beyond two-body perturbations [that is, studying the dynamics in the circular restricted three-body problem (CR3BP)], the Keplerian orbital elements are no longer well defined in certain regions of phase space, especially when the gravitational attractions of both the primary and secondary bodies have similar magnitudes, as occurs in the vicinity of the libration points. In this work, we define a generalization of orbit elements that can be applied in these regions and others. Specifically, we define a set of semi-analytical action-angle orbital elements that are defined locally about any bounded special solution in the CR3BP: equilibria, periodic orbits, and quasi-periodic orbits. Local action-angle orbital elements are defined using action-angle coordinates in the Birkhoff–Gustavson normal form about the bounded invariant manifold. We include detailed examples around the five equilibria $L_{1 - 5}$ in the Earth–Moon CR3BP.

I. Introduction

Orbital elements help organize a complicated dynamical environment by splitting it up into its fundamental parts. In the two-body problem, the Keplerian orbital elements are parameters that uniquely describe a trajectory: eccentricity, semimajor axis, inclination, right ascension of the ascending node, argument of periapsis, and true anomaly (or mean anomaly) [1]. The Keplerian orbital elements are used to describe the integrals of motion of the dynamical system [2]. Importantly, in unperturbed two-body dynamics, so-called Keplerian orbits are conic sections: circle, ellipse, parabola, and hyperbola. Once perturbations are imposed on two-body dynamics, the Keplerian orbital elements are no longer integrals of motion, and so one often studies the time evolution of the orbital elements or a sequence of osculating elements using the Lagrange planetary equations or Gauss equations [3]. Alternative orbital elements, such as the canonical Delaunay elements, exist [4]. Under additional perturbations to two-body dynamics, one may consider epicyclic orbital elements [5].

Moving in complexity beyond two-body perturbations [i.e., studying the dynamics in the circular restricted three-body problem (CR3BP)], the Keplerian orbital elements are no longer well defined in certain regions of phase space, especially when the gravitational attractions of both the primary and secondary bodies have similar magnitudes, as occurs in the vicinity of the libration points [6]. Whereas the relative two-body problem is integrable, the CR3BP is nonintegrable. So, unlike the conic sections of the two-body problem, the CR3BP has no global analytical solutions. Instead, the CR3BP emits special solutions: equilibrium points, periodic orbits, quasi-periodic orbits, and hyperbolic invariant manifolds (i.e., stable and unstable manifolds) [2]. These solutions inform our understanding of the dynamics of the CR3BP. Beyond the solutions themselves, the (linear) stability of the bounded special solutions (equilibria, and periodic and quasi-periodic orbits), either center or saddle, inform our understanding of the dynamics in the vicinity of the bounded special solutions. For example, the linearized dynamics about an equilibrium point in the CR3BP are integrable, and the frequencies corresponding to the center modes can be paired with a shooting method to compute nearby orbits. However, linear information in the vicinity of bounded special solutions has limited utility in explicitly and accurately describing the dynamics in such regions due to the low order of the expansion [7].

A normal form is a simplified form of the dynamics that is often used to elucidate the structure of solutions in a particular region of phase space [8–10]. For example, the linearized dynamics about an equilibrium point in the CR3BP can be put into a Birkhoff–Gustavson normal form via diagonalization [8,11]. Whereas the linear information cannot explicitly describe the dynamics in the vicinity of bounded special solutions, a higher-order expansion of the dynamics put into normal form can. Although it is known that integrals of time-independent Hamiltonian systems, in general (except $H$ ), do not exist [12], Gustavson found that by normalizing to an infinite order, additional formal integrals can be obtained [13]. In this so-called “indirect method,” one can apply the Lie series method, which is a semi-analytical technique, to bring the Hamiltonian into a Birkhoff–Gustavson normal form such that the integrals are immediately found [8,13–15]. The Lie series method is a semi-analytical technique used to construct a normal form, and hence compute approximate formal integrals while maintaining accuracy [16]. Using the Lie series method, one can construct a normal form and approximate formal integrals about any type of bounded special solution [17–19]. Jorba developed a computer algebra system to compute normal forms and approximate formal integrals of Hamiltonian dynamical systems, which is available in the public domain [17]; moreover, the computational methods of the work presented are based on Jorba’s computer algebra system, with modifications for the collinear libration points following Ref. [20] and extending these works to the spatial CR3BP. The primary contributions are the application and expansion of existing theory to apply the presented framework to orbit characterization in physical systems modeled using the CR3BP: in particular, in cislunar space.

An increase in interest of cislunar and translunar space has led to a need to describe the complex dynamics present in cislunar and translunar environments [21 –26]. Primarily, the lunar Gateway is planned to launch in 2024 and will enable a variety of robotic and manned missions to the Moon and beyond [22,24]. In addition, there are many missions past and present that have targeted or will target orbits in translunar space [25,26]. For example, missions such as Artemis (2007) and Genesis transited between $L_{1}$ and $L_{2}$ regions (Sun–Earth in the case of Genesis) [27–29]. Due to the presence of strong gravitational attraction from both the Earth and Moon in both cislunar and translunar space, mission analysts cannot readily rely on the typical two-body Keplerian orbital elements to describe the motion of a spacecraft in these regimes; moreover, using a more accurate model such as the CR3BP is required. The benefit of the Keplerian orbital elements is that they explicitly describe the motion of a spacecraft over all time. However, in the CR3BP, the dynamics are chaotic and trajectories are not as simply described as in the two-body setting. In this work, we develop the theory and a computational tool to describe CR3BP dynamics in cislunar and translunar space using local action-angle orbital elements, or action-angle elements for short, based on approximate integrals present near the equilibria of the CR3BP. For completeness, we extend the work to the remaining equilibria $L_{3 - 5}$ in the CR3BP.

The construction of these semi-analytical orbital elements in the CR3BP relies on a Birkhoff–Gustavson normal form about a bounded special solution: in particular, the approximate formal integrals of the normal form described by Gustavson [13]. Because this procedure is a local one, the so-called action-angle elements are locally defined. Action-angle elements rely on centering the coordinate system about the special solution, and so the treatment of equilibrium points, periodic orbits, and quasi-periodic orbits will be distinct in this initial step. For simplicity and clarity, we describe the construction about an equilibrium point. Action-angle elements explicitly and accurately describe the dynamics in the vicinity of a bounded special solution. In particular, the construction of action-angle elements defines regions of integrability in a nonintegrable dynamical system; in other words, half of the action-angle elements are approximate formal integrals of motion in the regions described. The other elements are corresponding angle coordinates (hyperbolic angle in the saddle case). Because the two-body orbital elements remain widely used in astrodynamics, the action-angle elements form a semi-analytical analog in the CR3BP and retain many of the same applications using a higher-fidelity model in the CR3BP. Furthermore, just as there is a map between Cartesian coordinates and Keplerian elements, there is a map between Cartesian coordinates and local action-angle orbital elements. We consider them to be “orbital elements” because they capture the motion of a body in a geometric way and are topologically similar to the two-body Delaunay action-angle variables [4].

Local parametrization of orbits in the CR3BP has been attempted before [21,30–32], but the tool presented is able to achieve a more intuitive parametrization to a higher order [33], considers a larger region than previous attempts, and incorporates bounded and unbounded special solutions [34]: namely, periodic orbits, quasi-periodic orbits, and hyperbolic invariant manifolds. In addition, this parametrization allows us a novel exploitation of the entire saddle manifold, beyond the typical stable and unstable manifolds, which are submanifolds of the saddle. We introduce this technique as a set of local action-angle orbital elements that can be used in a variety of mission analysis settings: some of which are discussed in the following.

The paper will be organized in the following way. First, we give a description of the CR3BP as a chaotic Hamiltonian dynamical system, including the series expansion of the Hamiltonian. Next, we construct the local action-angle orbital elements step by step (including discussions on a local Anosov splitting [35], the Birkhoff–Gustavson normal form, and action-angle coordinates for the center and saddle manifolds), and finally the action-angle elements themselves. We include an error analysis using several error metrics, including the radius of convergence, comparison to numerical integration, and isoenergetic level sets. We include detailed examples around the five equilibria $L_{1 - 5}$ in the Earth–Moon CR3BP divided into two groups based on their stability: the triangular points $L_{4, 5}$ , and the collinear points $L_{1 - 3}$ . In these sections, we describe the subtle differences in the action-angle elements, and we include several computations to show how action-angle elements explicitly describe the dynamics of the CR3BP in these regions. Additionally, we give explicit descriptions of particular mission applications for Genesis and Artemis-P1. Finally, we consider the implications of the work presented.

II. Dynamical Model: The Circular Restricted Three-Body Problem

In this work, we construct local canonical coordinates directly related to the trajectories in the CR3BP. To do this, we employ techniques from normal form theory that necessitate writing the Hamiltonian function of the system as

H = \sum_{n \geq 2} H_{n}

where

H_{n}

is a homogeneous polynomial of degree

n

. In this section, we describe the process of expanding the Hamiltonian of the CR3BP as an infinite series, as well as recentering the system about an equilibrium point. We also include a brief discussion on the dimensions of special solutions and how recentering the system depends on this discrepancy between types of bounded special solutions.

The CR3BP is a dynamical system in which we consider the motion of an infinitesimal particle under the gravitational forces of two massive primaries (modeled as point masses): in this case, the Earth and the Moon. We neglect the gravitational attraction of the infinitesimal particle on the primaries, and we assume the primaries are in a circular orbit about their mutual barycenter. The CR3BP is widely used in astrodynamics, and it can be used to study the motion of an asteroid under the attraction of Jupiter and the Sun or the motion of an artificial satellite in the Earth–Moon system, among many other systems [11].

We normalize the system by taking units of mass, length, and time such that the sum of masses of the primaries, the gravitational constant, and the period of the motion of the primaries is 1, 1, and $2 π$ , respectively [11]. We take $μ$ to be the normalized mass of the smaller primary, and $1 - μ$ is the normalized mass of the larger primary body so that $μ \in (0, 1 / 2]$ . We take the usual synodic reference frame: the origin is the mutual barycenter of the primaries, the $x$ axis is the line connecting the larger primary to the smaller primary, the $z$ axis is in the direction of the angular momentum of the primaries, and the $y$ axis completes the right-handed frame. Because motion takes place in this synodic frame, this system is autonomous by construction; i.e., it does not depend on time. Hence, defining momenta as $p_{x} = \dot{x} - y$ , $p_{y} = \dot{y} + x$ , and $p_{z} = \dot{z}$ , the equations of motion can be written in Hamiltonian form, with the Hamiltonian function given by [36]

H = \frac{1}{2} (p_{x}^{2} + p_{y}^{2} + p_{z}^{2}) + y p_{x} - x p_{y} - \frac{1 - μ}{r_{1}} - \frac{μ}{r_{2}}

(1)

where

r_{1}

is the distance from the larger primary, and

r_{2}

is the distance from the smaller primary. In summary, the CR3BP is an autonomous Hamiltonian system that can be used to model the motion of an artificial satellite in the Earth–Moon system.

Unlike the two-body problem (2BP), the CR3BP is not an integrable dynamical system; it is not generically solvable. Although the solutions of the 2BP are conic sections [1], there is no global integrable solution of the CR3BP, and trajectories can exhibit chaotic behavior [1,2]; however, special solutions exist: equilibria, periodic and quasi-periodic orbits, as well as hyperbolic invariant manifolds. Figure 1 shows the five equilibrium points in the Earth–Moon CR3BP. $L_{1 - 3}$ are called the collinear equilibria, and $L_{4 - 5}$ are called the triangular equilibria. We compute the location of the triangular equilibria explicitly and the collinear points from the Euler quintic equation [36], given by distance $γ_{j}$ for $j = 1, 2, 3$ . There exist periodic and quasi-periodic orbits emanating from these equilibria. Associated to unstable bounded special solutions, there are hyperbolic invariant manifolds. In this work, we show how the action-angle elements exploit the whole saddle manifold rather than the usual submanifolds of the saddle, i.e., the stable and unstable manifolds.

Local action-angle orbital elements are defined locally about any bounded special solution in the CR3BP: equilibria, periodic orbits, and quasi-periodic orbits. To define this set of local coordinates, we must center our coordinates about the particular special solution. We treat the cases of the collinear equilibria and triangular equilibria separately to avoid abundant notation.

Fig. 1
The five equilibrium points of the Earth–Moon CR3BP.

A. Recentering About Equilibria

In the case of the triangular equilibria $L_{4, 5}$ , the Hamiltonian of the CR3BP centered at the equilibria is given by [17]

H = \frac{1}{2} (p_{x}^{2} + p_{y}^{2} + p_{z}^{2}) + y p_{x} - x p_{y} + (\frac{1}{2} - μ) x \pm \frac{\sqrt{3}}{2} y - \frac{1 - μ}{r_{1}} - \frac{μ}{r_{2}}

(2)

where

r_{1}

and

r_{2}

are defined as shown earlier in this paper, “

+

” is taken for

L_{4}

, and “

-

” is taken for

L_{5}

. Because the linear stability of the coordinate origin affects the forms of the action-angle elements, the mass ratio is taken below the Routh critical value (

μ \leq μ_{R} \approx 0.0385

) so that the origin is linearly stable [1].

In the case of the collinear equilibria, this is the result of the following translation and scaling:

X = \mp γ_{j} x + μ + a, Y = \mp γ_{j} y, Z = γ_{j} z

(3)

where the upper sign corresponds to

L_{1, 2}

, the lower one to

L_{3}

a = - 1 + γ_{1}

for

L_{1}

a = - 1 - γ_{2}

for

L_{2}

, and

a = γ_{3}

for

L_{3}

. In the case of periodic or quasi-periodic orbits, this translation is more nuanced and forces the system to be nonautonomous [37], i.e., time dependent. This case is worthy of deeper exploration and will be the topic of future work.

B. Series Expansion of Hamiltonian

Next, in order to perform the necessary Lie transformations to put the Hamiltonian function into its Birkhoff–Gustavson normal form at a particular order, we must write

H = \sum_{n \geq 2} H_{n}

where

H_{n}

is a homogeneous polynomial of degree

n

. The only terms of

H

not already in this form are the dipole terms:

1 - μ / r_{1}

and

μ / r_{2}

. We use the following dipole expansion to expand these terms:

\frac{1}{\sqrt{(x - A)^{2} + (y - B)^{2} + (z - C)^{2}}} = \frac{1}{D} \sum_{n = 0}^{\infty} (\frac{ρ}{D})^{n} P_{n} (\frac{A x + B y + C z}{D ρ})

(4)

where

D^{2} = A^{2} + B^{2} + C^{2}

ρ^{2} = x^{2} + y^{2} + z^{2}

, and

P_{n}

is the Legendre polynomial of degree

n

. Using the preceding expansion, we obtain the original Hamiltonian now recentered about an equilibrium point and written as an infinite series:

H = \frac{1}{2} (p_{x}^{2} + p_{y}^{2} + p_{z}^{2}) + y p_{x} - x p_{y} - \sum_{n \geq 2} c_{n} (μ) ρ^{n} P_{n} (\frac{x}{ρ})

(5)

where

c_{n} (μ)

are real coefficients depending on

μ

and

γ_{j}

given by Jorba and Masdemont [11].

III. Local Action-Angle Orbital Elements

The main contribution of this work is the construction of local action-angle orbital elements in the CR3BP. In this section, we construct the action-angle elements by using an equilibrium point as the center of our coordinate system for clarity. We begin with the preliminaries for Birkhoff–Gustavson normal form construction: that is, linearization, diagonalization, and complexification. We then discuss in brief the Lie series method, citing relevant references for further investigation by the interested reader. Next, we describe the final coordinate transformation to send normal form coordinates into action-angle elements. In light of the main contribution of this work, we provide a high-level overview of the method for clarity of future implementation. Finally, we include remarks about extending this procedure to higher-dimensional bounded special solutions.

A. Preliminaries

In the previous section, we described the initial recentering of the coordinate system about an equilibrium point. To construct a higher-order Birkhoff–Gustavson normal form about an equilibrium point, we first put the linear terms into their complex diagonal form. This procedure explicitly shows the integrability of the linear dynamics, and consequently illustrates the linear stability. Furthermore, this form makes the normal form computation more efficient [17].

First, we linearize the dynamics by taking only the second-order terms in the series expansion of Eqs. (1) and (2). Next, one can construct a real linear symplectic change of variables such that the Hamiltonian [either Eqs. (1) or (2)] is reduced to its (real) normal form [17]. This computation for both the collinear and triangular equilibria can be found in the literature: most notably in the appendices of [17]. In this way, we put the quadratic terms (which correspond to the linear dynamics) into real normal form:

H_{2}^{L_{1 - 3}} = λ x p_{x} + \frac{ω_{1}}{2} (y^{2} + p_{y}^{2}) + \frac{ω_{2}}{2} (z^{2} + p_{z}^{2})

(6)

H_{2}^{L_{4 - 5}} = \frac{ω_{1}}{2} (x^{2} + p_{x}^{2}) + \frac{ω_{2}}{2} (y^{2} + p_{y}^{2}) + \frac{ω_{3}}{2} (z^{2} + p_{z}^{2})

(7)

where the equations correspond to the collinear and triangular points, respectively; and the coefficients

λ

ω_{j} \in R

for

j = 1, 2, 3

are constant coefficients depending on the equilibrium point of interest

L_{1 - 5}

and the mass parameter

μ

. We see that

H_{2}

differs, depending on the linear stability.

Throughout this work, we use the complex normal form for $H_{2}$ because it is convenient for constructing the normal form, as we will see later. This complexifying change is dependent on the linear stability of the equilibrium point. We need only complexify the center parts of the real linear normal form because the saddle term is already in the appropriate form. The complexification of the center terms is given by

x = \frac{q_{1} + i p_{1}}{\sqrt{2}}, p_{x} = \frac{i q_{1} + p_{1}}{\sqrt{2}}

(8)

where (

x

p_{x}

) can be replaced with (

y

p_{y}

) or (

z

p_{z}

), depending on the application, so long as the right-hand side is changed to (

q_{2}

p_{2}

) or (

q_{3}

p_{3}

), respectively. In this way, Eqs. (1) and (2) have been reduced to their linear complex normal forms, respectively:

H_{2}^{L_{1 - 3}} = λ q_{1} p_{1} + i ω_{1} q_{2} p_{2} + i ω_{2} q_{3} p_{3}

(9)

H_{2}^{L_{4 - 5}} = i ω_{1} q_{1} p_{1} + i ω_{2} q_{2} p_{2} + i ω_{3} q_{3} p_{3}

(10)

where

λ

ω_{j} \in R

for

j = 1, 2, 3

are all positive real numbers.

Now, we could define the action-angle variables: $I_{j} = q_{j} p_{j}$ and $θ_{j}^{c} = \arctan (p_{j} / q_{j})$ or

θ_{j}^{s} = \ln \frac{\sqrt{q_{j}}}{\sqrt{p_{j}}}

for

j = 1, 2, 3

, where

θ_{j}^{c}

and

θ_{j}^{s}

correspond to the center and saddle, respectively. One can check that

I_{j}

are integrals of the linear normal form

H_{2}

, and

θ_{j}

are the symplectic conjugate variables, i.e., variables such that

(q, p) \mapsto (I, θ)

is symplectic; see the Appendix for more on the hyperbolic action-angle variables. Hence, the linear normal form is integrable with integrals of motion corresponding directly to the center and saddle manifolds. In these action-angle variables, the equations of motion become constant actions and linear angles. The goal of this work is to include higher-order terms in the Hamiltonian while maintaining integrability of the dynamics. As such, we will construct a higher-order Birkhoff–Gustavson normal form that preserves the linear integrals so that they become approximate formal integrals of the series expansion. When we compute the higher-order normal form, we will use the preceding action-angle coordinate transformations after computing the normal form rather than after the linearization, as noted earlier in this paper.

In the earlier part of this paper, the linear dynamics have been split into their center and saddle directions via composition of canonical maps: position momenta, diagonalization (from the appendices of Ref. [17]), and complexification. Because canonical transformations are diffeomorphisms, we have made a diffeomorphism of the configuration manifold $R^{3}$ and its tangent space to itself with clearly marked directions of expansion (i.e., unstable manifold), contraction (i.e., stable manifold), and nonexpanding and noncontracting. In mathematics, such a transformation is referred to as a local Anosov splitting of the phase space at each equilibrium point in the system. Describing the flow of the CR3BP as a local Anosov flow near bounded special solutions (i.e., equilibria in the paper presented) greatly clarifies the dynamics. In the following section, we will show how to expand the region of validity for the local Anosov approximation.

B. Lie Series Method

Next, we will give a brief description of the Birkhoff–Gustavson normal form. Because the primary contribution is not the computation of the normal form, we will refer the reader to fundamental works, such as that by Jorba [17], where appropriate. We will then describe step by step the construction of local action-angle orbital elements so that the procedure and implementation are as clear as possible. We must also include discussion on the differences in action-angle element constructions for each stability type, e.g., center and saddle. We refer to the Appendix for more details on hyperbolic action-angle coordinates, which are the action-angle elements corresponding to the saddle.

We now give a high-level overview of the Lie series method because it is used to compute a Birkhoff–Gustavson normal form about the equilibria in the CR3BP. For more on this technique, see Refs. [2,9,16,17,20,38]. The Lie series method is a technique used often to simplify the form of the dynamics canonically while maintaining accuracy. The Lie series method relies on the Lie transform, which is the flow for $t = 1$ along a vector field defined by a generating function. This procedure relies on the following two facts. First, if $Φ_{t} (q, p)$ is a time $t$ flow of a Hamiltonian system, then $(Q, P) = Φ_{t} (q, p)$ is a canonical transformation. Second, the time derivative of any smooth function $F (q, p)$ in a Hamiltonian system governed by Hamiltonian function $G (q, p)$ is the Poisson bracket of $F$ with $G$ :

\frac{d}{d t} F = {F, G}

where {

F

G

} is the Poisson bracket of functions defined by

{F, G} = \sum_{j = 1}^{n} (\frac{\partial F}{\partial q_{j}} \frac{\partial G}{\partial p_{j}} - \frac{\partial F}{\partial p_{j}} \frac{\partial G}{\partial q_{j}})

Note that if $F$ is a polynomial of degree $r$ and $G$ is a polynomial of order $s$ , then the Poisson bracket { $F$ , $G$ } is a polynomial of degree $r + s - 2$ .

Now, let $F (t) = H \circ Φ_{t} (q, p)$ , where $Φ_{t} (q, p)$ is the time $t$ flow of a Hamiltonian system $G (q, p)$ , and $H$ is the Hamiltonian of the CR3BP centered at an equilibrium point given in Eq. (5). Performing a Taylor series expansion of $F$ about $t = 0$ gives

F (t) = F (0) + F^{'} (0) t + \frac{1}{2!} F^{''} (0) t^{2} + \dots + \frac{1}{n!} F^{(n)} (0) t^{n} + \dots

(11)

Evaluating the preceding expansion at $t = 1$ gives

F (1) = H \circ Φ_{1} (q, p) = H + {H, G} + \frac{1}{2!} {{H, G}, G} + \dots + \frac{1}{n!} L_{G}^{n} (H) + \dots

(12)

where

L_{G}^{n} (H)

is the

n

th Lie derivative of

H

generated by

G

, i.e., the

n

-fold nested Poisson bracket. We denote this expression as

\hat{H} \equiv F (1)

. In short, what we have done is define a canonical transformation by defining a new vector field (generated by

G

) on the phase space and flowing along for

t = 1

. The function of

G

generating this new vector field is freely chosen for each application.

The Lie series method also provides a direct change of variables, following a similar form to Eq. (12):

{\hat{q}}_{j} = q_{j} + {q_{j}, G_{k}} + \frac{1}{2!} {{q_{j}, G_{k}}, G_{k}} + \dots + \frac{1}{n!} L_{G_{k}}^{n} (q_{j}) + \dots

(13)

{\hat{p}}_{j} = p_{j} + {p_{j}, G_{k}} + \frac{1}{2!} {{p_{j}, G_{k}}, G_{k}} + \dots + \frac{1}{n!} L_{G_{k}}^{n} (p_{j}) + \dots

(14)

for

j = 1, 2, 3

and

k = 3, \dots, N

. The inverse transformation using the Lie series method is geometrically intuitive. We obtain the direct transformation by flowing along the vector field generated by

G

for a time of

t = 1

_; to obtain the inverse transformation, we flow along the vector field generated by

- G

for a time of

t = 1

(or, equivalently, flow along the vector field generated by

G

for a time of

t = - 1

). That is, the inverse transformation is given by the preceding expressions by replacing

G_{k}

with

- G_{k}

Because we seek to construct a higher-order Birkhoff–Gustavson normal form, we will choose $G$ at each order such that $\hat{H}$ preserves the integrability of $H_{2}$ . To do this, it suffices to choose $G$ so that we eliminate all terms at each order. Intuitively, the integrable $H_{2}$ remains unchanged by each transformation, and we obtain a higher-order normal form through a composition of transformations. We will show in the following why this cannot happen in this case, as well as how we are able to work around this problem.

First, we consider eliminating third-order terms. By Eq. (12), the third-order terms in the transformation are

{\hat{H}}_{3} = H_{3} + {H_{2}, G_{3}}

(15)

where

H_{2}

and

H_{3}

are known [given in Eq. (5)],

{\hat{H}}_{3} = 0

because we seek to eliminate third-order terms in the transformed Hamiltonian, and we solve this homological equation for

G_{3}

. We can do this explicitly:

G_{3} (q, p) = \sum_{| k_{q} | + | k_{p} | = 3} \frac{- h_{k_{q}, k_{p}}}{i ⟨ k_{p} - k_{q}, ω ⟩} q^{k_{q}} p^{k_{p}}

(16)

where we use the multi-indices

k_{q}

and

k_{p}

, and

ω

is the linear frequency vector corresponding to the equilibrium point. Performing this transformation (i.e., moving along the flow generated by

G_{3}

for unit time) gives the following transformed Hamiltonian:

\hat{H} = H_{2} + 0 + {\hat{H}}_{4} + {\hat{H}}_{5} + \dots

(17)

Now, if we take the fourth-order terms and attempt the same procedure, we obtain

G_{4} (q, p) = \sum_{| k_{q} | + | k_{p} | = 4} \frac{- h_{k_{q}, k_{p}}}{i ⟨ k_{p} - k_{q}, ω ⟩} q^{k_{q}} p^{k_{p}}

(18)

Because $| k_{q} | + | k_{p} | = 4$ , or an even number, this generating function cannot exist because the denominator is zero when $k_{q} = k_{p}$ . So, on the one hand, we must keep those terms in ${\hat{H}}_{4}$ , and so on, for all even terms. On the other hand, we are always able to eliminate all odd-ordered terms.

Performing this procedure up to some order $N$ , we obtain a simplified Hamiltonian of the form

\hat{H} = H_{2} + {\hat{H}}_{N} (q_{1} p_{1}, q_{2} p_{2}, q_{3} p_{3}) + R_{N} (q, p)

(19)

where

{\hat{H}}_{N}

includes terms of order

N

and all even orders (except two) less than or equal to

N

, and where

R_{N}

is the remainder term, i.e., terms of order greater than

N

. Note that

H_{2}

and

{\hat{H}}_{N}

are functions of three variables,

q_{j} p_{j}

for

j = 1, 2, 3

. As in the linear normal form, defining

I_{j} = q_{j} p_{j}

yields three integrals of motion (in the complex variables). We transform the complex variables into real variables (

q

p

) using the inverse of Eq. (8) to obtain the real Hamiltonian:

\hat{H} = H_{2} + {\hat{H}}_{N} (q_{1} p_{1}, \frac{1}{2} (q_{2}^{2} + p_{2}^{2}), \frac{1}{2} (q_{3}^{2} + p_{3}^{2})) + R_{N} (q, p)

(20)

where we follow the preceding notation in the complex normal form. Observe that the forms of the integrals may change but that

{\hat{H}}_{N}

remains a function of only the integrals. In this way, we say that

\hat{H}

is integrable up to order

N

. The Hamiltonian function

\hat{H}

in Eq. (20) is called the real Birkhoff–Gustavson normal form of order

N

C. Local Action-Angle Orbital Elements

Once the real Birkhoff–Gustavson normal form is computed up to order $N$ , we define action-angle variables as in a previous discussion, depending on the linear stability of the equilibrium point. The center and saddle integrals and conjugate angles are defined by

I_{c} = \frac{1}{2} (q^{2} + p^{2}), θ_{c} = \arctan (p / q)

(21)

I_{s} = q p, θ_{s} = \ln \frac{\sqrt{q}}{\sqrt{p}}

(22)

respectively. For

I_{s} > ε

, the action-angle elements are given by action-angle variables corresponding to the linear stability. We consider the cases of triangular and collinear equilibria separately, and we consider the case of

I_{s} > ε

for sufficiently small

ε > 0

(see the Appendix for

I_{s} < ε

For the triangular equilibria, assuming a mass parameter of $μ \leq μ_{R} \approx 0.0385$ (where $μ_{R}$ is Routh’s critical value [1]), the linear stability is center $\times$ center $\times$ center. Hence, the action-angle elements in these regions are given by ( $I_{1}$ , $θ_{1}$ , $I_{2}$ , $θ_{2}$ , $I_{3}$ , $θ_{3}$ ), where ( $I_{j}$ ) for $j = 1, 2, 3$ are computed using ( $I_{c}$ , $θ_{c}$ ), respectively, using the order $N$ normal form coordinates, $q_{j}$ and $p_{j}$ .

For the collinear equilibria, the linear stability is saddle $\times$ center $\times$ center. Hence, for $I_{s} > ε$ , the action-angle elements in these regions are given by ( $I_{1}$ , $θ_{1}$ , $I_{2}$ , $θ_{2}$ , $I_{3}$ , $θ_{3}$ ), where ( $I_{1}$ , $θ_{1}$ ) are computed using ( $I_{s}$ , $θ_{s}$ ), and ( $I_{j}$ , $θ_{j}$ ) for $j = 2, 3$ are computed using ( $I_{c}$ , $θ_{c}$ ) from earlier in this paper by using the order $N$ normal form coordinates, $q_{j}$ and $p_{j}$ . Now, the case of $I_{s} < ε$ is discussed in greater detail in the Appendix. In short, for sufficiently small $I_{s}$ , we neglect the action-angle transformation and retain the standard normal form coordinates ( $q_{1}$ , $p_{1}$ ) corresponding to the saddle.

Ultimately, by composing the coordinate transformations above, we have obtained a sequence of coordinate transformations that can be succinctly represented as

(x, y, z, \dot{x}, \dot{y}, \dot{z}) \mapsto (I_{1}, θ_{1}, I_{2}, θ_{2}, I_{3}, θ_{3})

(23)

defined in the regions around each equilibrium point in the CR3BP. Note that, because scaling does not preserve volume, the preceding transformation is not canonical, although it is canonical once scaled and recentered about the equilibrium point. Furthermore, note that the inverse transformation taking action-angle elements to Cartesian coordinates also exists within the region. Because the Lie series method is ultimately a composition of asymptotic series expansions, it is a divergent procedure [39]; therefore, there will be an optimal finite order to take the Birkhoff–Gustavson normal form.

D. High-Level Overview

The purpose of this paper is to introduce semi-analytical orbital elements for multibody systems. For this purpose, we include here a high-level overview of the computational procedure for application of action-angle elements about an equilibrium point, leaving the application periodically perturbed systems as a topic for future work. Note that each step varies in degree of difficulty, in implementation as well as computational time:

Transform $(x, y, z, \dot{x}, \dot{y}, \dot{z}) \mapsto (x, p_{x}, y, p_{y}, z, p_{z})$ .
Scale and recenter about equilibrium point of interest.
Diagonalize using a linear symplectic change of variables.
Complexify center terms.
Perform successive Lie transformations and realify variables.
Perform action-angle coordinate transformations, depending on the linear stability of the equilibrium point.

IV. Error Analysis

The methods applied in this work are semi-analytical normal form techniques dealing with series truncation and successive coordinate transformations, and so some error analysis is warranted. The acceptable level of error is dependent on the application: a computer-assisted proof would require greater precision than a simple heuristic survey of regional dynamics in the CR3BP. In this section, we discuss in brief several metrics for error analysis: radius of convergence, comparison with numerical integration over time, and divergence in energy level sets. We show that radius of convergence is the most strict error metric of the three, whereas comparison via numerical integration is useful for individual trajectories and divergence in energy level sets exploits the reduced dimension of the integral space to determine maximal bounds of validity. We begin with some definitions.

Following Jorba et al. [17,37], in order to have an idea of the radius of convergence of the series in Birkhoff–Gustavson normal form, we have computed the values

r_{n} = \sqrt[n]{‖ H_{n} ‖_{1}}, where ‖ H_{n} ‖_{1} = \sum_{| k | = n} | h_{k} |, 3 \leq n \leq N

(24)

with

h_{k}

as the coefficient of the monomial of exponent

k

. Using the root test given earlier in this paper, we take

1 / r

normalized such that if

1 / r = 1

, the radius of convergence is precisely the distance from the libration point to the nearest singularity (Moon for

L_{1, 2}

, and Earth for

L_{3}

). Figure 2 shows the normalized ratio and root tests up to order 32 for

L_{1 - 3}

. Note the tradeoff in convergence the radius and order of approximation.

Although the radius of convergence is a useful tool for pure mathematical applications, this method generally underestimates the region of approximation for astrodynamics applications using a high-order expansion. Conversely, although the radius of convergence may be very high for a low-order expansion, it is not necessarily accurate out to its radius of convergence. In other words, as the order of the series expansion increases, the radius of convergences decreases, although we may expect the series to increase in accuracy as more terms are included. So, we consider additional metrics, focusing on the $L_{1}$ point for the remainder of the analysis.

Fig. 2
Normalized radius of convergence for the collinear equilibria.

The simplest such metric is a direct comparison between normal form integration and numerical integration (via the Runge–Kutta method) of the exact differential equations. We compare the integration over several amplitudes for a fixed time step. We take initial conditions $I_{2} = λ_{0}$ and $I_{1, 3} = 0$ with initial phases $θ_{1 - 3} = 0$ for $t = 0$ , save the corresponding solution at $t_{f} = 0.01$ , and transform both points to synodical coordinates to obtain two points: $v_{0}$ and $v_{1}$ . Then, we integrate via the Runge–Kutta method in the CR3BP from $v_{0}$ until $t_{f}$ , and we denote this point $v_{0}^{1}$ . Table 1 gives the difference of $‖ v_{1} - v_{0}^{1} ‖_{2}$ . Note that the local error of the numerical integration is of the order of $10^{- 16}$ , and the normal form (and the corresponding changes of variables) has been computed up to degree 32.

We can also use numerical integration to compare the dynamical effect on a single trajectory over time. We consider a planar Lyapunov orbit about $L_{1}$ with $I_{2} = 0.5$ and $I_{1, 3} = 0$ . Starting from $θ_{2} = 0$ , we integrate for 10 nondimensional time units (approximately 43.4 days). Figure 3 shows the trajectories computed using the normal form and Runge–Kutta integration, with the difference shown in log scale in the following. For Figs. 3 –5, we denote by NF the normal form integration, and by RTBP the Runge–Kutta integration of the restricted three-body problem equations of motion. Note that the error grows over time and diverges after several periods. Considering the transformation of the CR3BP numerically integrated trajectory into action-angle elements, we observe that the instability of the trajectory is captured by the unstable manifold component; see Figs. 4 and 5.

One can additionally compare the action-angle elements of a numerically computed (for example, via shooting method) periodic orbit. In the following, we include two cases of this. The first case is a planar Lyapunov periodic orbit. Because this is within our approximation region, we expect the action-angle elements to have an approximately zero saddle integral, a constant planar integral, and a linearly changing planar angle. We observe this in Fig. 6. In the second case, we consider a northern halo orbit, shown in Fig. 7; notably, this orbit falls outside of our approximation region because we have killed the resonant terms of the Hamiltonian expansion that account for the bifurcation that gives rise to the halo orbit family [40]. Hence, there is uncertainty as to the action-angle element representation for this orbit. Figure 7 shows the action-angle elements over one orbit period. The saddle integral remains bounded below $10^{- 7}$ , with no sign of stable/unstable drift; the planar and vertical integrals, $I_{2}$ and $I_{3}$ , vary on the order of $10^{- 3}$ in a symmetric way; the planar and vertical angles, $θ_{2}$ and $θ_{3}$ , are approximately linear in time. This example indicates there are trajectories theoretically not accounted for by the Hamiltonian expansion, approximation, and truncation; but they maintain a local action-angle orbital element representation, albeit at a lower precision. Figures 4 and 8 are reminiscent of the semimajor axis and mean motion of two-body Keplerian orbits over time under additional perturbations.

Fig. 3
Comparing numerical integration and normal form integration at order 32 about $L_{1}$ .

Fig. 4
Action-angle elements for a planar orbit and differences between integrated trajectories.

Fig. 5
Saddle manifold elements (left) and stable/unstable manifold normal form coordinates (right)

The downside of comparing via integration is that it does not give a radial estimate defining a region of approximation to use. As a preliminary measure, to understand this region for bounded orbits, one can compute isoenergetic level sets of $H$ to see where divergence begins to occur, as shown in Fig. 9. The lobes are due to the series truncation and form the boundary; as the order of the series increases, the lobes shrink to the origin, as seen in figures 5.12 and 5.13 in the work of Giorgilli [39]; however, approximately linear behavior can be considered useful, depending on the application. For example, the $h = 0.7$ level set could be used in the region of $‖ I ‖ \leq 0.25$ (i.e., the ball with a radius of 0.25 in the $I_{2} - I_{3}$ plane) because this avoids the disconnected lobes formed near the edges. Moreover, as Fig. 9 illustrates, in the lower-dimensional integral space, we can identify points in the $I_{2} - I_{3}$ plane with periodic and quasi-periodic orbits. Although not a rigorous result, this method uses a major benefit of the Birkhoff–Gustavson normal form: visualizing the six-dimensional phase space in the two- or three-dimensional integral space.

Fig. 6
Local action-angle orbital elements for a numerically computed $L_{1}$ planar Lyapunov orbit.

Fig. 7
Numerically computed Northern halo orbit about $L_{1}$ .

Fig. 8
Local action-angle orbital elements for a numerically computed $L_{1}$ Northern halo orbit.

Furthermore, extending Fig. 9 to include a nonzero saddle integral shows the isoenergetic level surfaces in the full three-dimensional integral space. Here, points correspond to saddle surfaces emanating off of quasi-periodic orbits. Because these are generically volumes in the configuration space, we show in Fig. 10 an example of a saddle trajectory corresponding to a point in the integral space, imposing the condition that $θ_{j} = 0$ for $j = 2, 3$ , thus picking out a single trajectory. Note that the trajectory is unbounded, although it comes close to the Lissajous orbit, because we have chosen the saddle integral to be small.

Although analyzing the error of this approach is necessary, the coordinate transformations between synodical coordinates and action-angle elements exist and can be used interchangeably to identify arbitrary points in the phase space with their bounded orbit and saddle trajectory counterparts.

Fig. 10
Iso-energetic level surfaces of $H$ about $L_{1}$ .

Note that the expansion order presented is neither necessarily nor likely the optimal order of expansion. Depending on the application, the optimal order may be lowered; moreover, because lower-order transformations can be recovered from any higher-order computation, the high order presented is used to show the computational capabilities of the method.

V. Application to Triangular Equilibria

In this and the following section, we include numerical examples of local action-angle orbital elements computed in the vicinity of the five equilibria in the Earth–Moon CR3BP. In this section, we discuss the triangular equilibria, $L_{4}$ and $L_{5}$ , which have linear stability of the following type: center $\times$ center $\times$ center. This means that in the vicinity of the equilibria, there is a six-dimensional center manifold. In a previous section, we described the construction of action-angle elements for this stability case; here, we study two- and three-dimensional quasi-periodic invariant tori in these regions.

Figure 11 shows an example computation illustrating a two-dimensional planar quasi-periodic orbit about $L_{4}$ . For this trajectory, the values of integrals $I_{1}$ and $I_{2}$ are fixed at the same nonzero value of $I_{1} = I_{2} = 10^{- 4}$ and $I_{3} = 0$ , and so there is no out-of-plane oscillation. The initial angles are set to zero, and the trajectory is integrated numerically using the Birkhoff–Gustavson normal form. Note that the angle trajectory jumps back to $θ_{j} = 0$ whenever it crosses $θ_{j} = 2 π$ . This figure was used to verify the numerical procedure with the results of Jorba [17]. The left-side image shows the trajectory in the configuration space, and the right-side image shows the trajectory in action-angle elements space: in other words, the angle space (because the actions remain constant). We see that the figure on the right exhibits linear motion, which is far simpler than the trajectory in the figure on the left.

Using the local action-angle orbital elements, the integrals (actions) split the phase space into its fundamental modes of motion. In the case of $L_{4 - 5}$ , there are three center integrals. These integrals, along with their symplectic conjugate angle coordinates, can be used to better understand high-dimensional dynamical solutions of the system. For example, three-dimensional quasi-periodic invariant tori exist in several regions of the Earth–Moon system; however, these orbits have long been difficult to visualize because they fill an entire volume in the configuration space. Figure 12 shows an example computation of a three-dimensional spatial quasi-periodic orbit about $L_{4}$ . Similar to the setup of Fig. 11, we numerically integrate via Birkhoff–Gustavson normal form the elements $I_{j} = 10^{- 4}$ for $j = 1, 2, 3$ , and the initial angles are all set to zero. The resulting trajectory appears as a bounded chaotic mess in the configuration space (left). Yet, in the angle space, the trajectory is simply linear (right). Hence, the local action-angle orbital elements can be used to simplify high-dimensional dynamical solutions of the system.

Fig. 11
Two-dimensional planar quasi-periodic orbit about Earth–Moon $L_{4}$ .

Fig. 12
Three-dimensional quasi-periodic orbit about Earth–Moon $L_{4}$ .

VI. Application to Collinear Equilibria

Similar to the previous section, we include numerical examples about the collinear equilibria $L_{1 - 3}$ , which all have linear stability of saddle $\times$ center $\times$ center. This means that there is a four-dimensional center manifold and a two-dimensional saddle manifold. The saddle manifold can be parametrized by the combination of the stable and unstable manifolds. In a previous section, we constructed action-angle elements for this stability case and described the modifications from the triangular equilibria for this construction. We use these elements to study the dynamics near $L_{1 - 3}$ .

First, we can validate this numerical procedure by computing invariant objects about each of the collinear equilibria. In particular, Fig. 13 shows a section of the family of planar Lyapunov periodic orbits about $L_{3}$ . These orbits are computed by generating a list of actions of $I_{2} \in [0, 0.05]$ and for each action generating a list of points along the orbit of $θ_{2} \in [0, 2 π)$ . We compute out approximately to the boundary of our valid region of approximation.

Next, Fig. 14 shows the stable and unstable manifolds of a planar Lyapunov periodic orbit about $L_{1}$ . These manifolds are numerically integrated using the full dynamics of the CR3BP with a Runge–Kutta method, stopping once the $x$ value reaches the Earth or the Moon, and the manifolds are seeded by perturbing each point of the periodic orbit along the stable or unstable manifold. This is a familiar procedure because we often perturb along the stable or unstable eigenvector; however, using action-angle elements, we have higher-order approximations for the stable and unstable manifolds, and the perturbation is built into the coordinates ( $q_{1}$ , $p_{1}$ ) corresponding to the unstable and stable manifolds, respectively. For example, fix $I_{2} = I_{2}^{*}$ and let $θ_{2} \in [0, 2 π)$ . Along the periodic orbit, the remaining coordinates ( $I_{1}$ , $θ_{1}$ , $I_{3}$ , $θ_{3}$ ) or ( $q_{1}$ , $p_{1}$ , $I_{3}$ , $θ_{3}$ ) are all zero. To seed the stable and unstable manifold points then, we simply set $p_{1}$ or $q_{1}$ to some $ε > 0$ along each point on the periodic orbit (parametrized by $θ_{2}$ ).

A persistent problem in preliminary mission design, especially with orbits as abundant as quasi-periodic orbits, includes categorization and classification. In particular, two-dimensional (2-D) quasi-periodic orbits in the CR3BP lie in two-parameter families. The increase in dimension from periodic orbits to quasi-periodic orbits makes the latter more challenging to catalog as compared to the former. On one hand, to catalog spatial periodic orbits, one needs only seven pieces of data: a six-dimensional initial state in Cartesian coordinates, and a period (or Jacobi constant). On the other hand, to catalog quasi-periodic orbits, the information required is dependent on the method used to compute invariant tori, and more data are required, regardless of the method. For instance, to catalog a 2-D quasi-periodic orbit computed by using a flow map method such as the one developed by Gomez and Mondelo, later modified by Olikara and Scheeres [41], the following data are required: multiple six-dimensional states along an invariant curve (often, multiple curves along the orbit, the stroboscopic time, and the rotation number).

Fig. 14
Stable/unstable manifolds of a periodic orbit seeded via Birkhoff–Gustavson normal form.

The use of local action-angle orbital elements can be a complementary method of cataloging quasi-periodic orbits because they give a local parametrization of families of orbits within a given region. Using action-angle elements, a mission design engineer is able to classify orbits quickly and with ease due to the simple form of the series expansion. Figure 15 shows a sample of $L_{2}$ Lissajous orbits with varying planar oscillation of $I_{2} \in [0, 0.2]$ and fixed vertical oscillation of $I_{3} = 0.1$ . Note that one can understand Fig. 15 via the isoenergetic level sets in Fig. 9: we fix $I_{3} = 0.1$ and move along the $I_{2}$ axis to generate Lissajous tori. Furthermore, by using action-angle elements, we exploit coordinates split into their fundamental dynamical modes to a high degree in order to classify the two-parameter family of $L_{2}$ quasi-periodic Lissajous orbits.

Although Figs. 13 –15 are more familiar images, albeit computed by different means, they show either the fixed bounded invariant objects (periodic and quasi-periodic orbits) or hyperbolic invariant manifolds (stable and unstable manifolds) of a particular orbit. Now, there is a higher-dimension version of the hyperbolic invariant manifolds, namely, the saddle manifold. Whereas the stable and unstable manifolds add one dimension to that of their bounded orbit counterparts (e.g., periodic orbit or quasi-periodic orbit), the saddle manifold adds two dimensions, with the stable and unstable manifolds locally spanning the saddle manifold. One of the contributions of local action-angle orbital elements is giving dynamically meaningful coordinates to arbitrary points in phase space near a bounded special solution. In this work, these regions are limited to the vicinity of the equilibria; even so, we gain insight into the dynamics of the Earth–Moon system and transitions between regions.

Fig. 15
Collection of Earth–Moon $L_{2}$ Lissajous quasi-periodic orbits at a fixed vertical action.

Every point that passes through an appreciable neighborhood of any bounded special solution in the CR3BP can be represented in its fundamental dynamical modes. For example, using action-angle elements, we can study the saddle manifold locally in the following way. The saddle manifold is locally diffeomorphic to the $q_{1} - p_{1}$ plane, with the stable and unstable manifolds on the $p_{1}$ and $q_{1}$ axes, respectively. Trajectories that do not lie on these axes are represented using the saddle integral $I_{1}$ and the symplectic conjugate angle $θ_{1}$ , where $θ_{1} \in (- \infty, \infty)$ and $θ_{1} = 0$ represent the points of closest approach. The quadrant is determined by the imaginary offset of $θ_{1}$ _; see the Appendix for more details. If this offset is an integer multiple of $π$ , then the point is in either the first or third quadrant of the saddle manifold. In these quadrants, if $θ_{1} = 0$ , then the point in phase space is crossing the so-called “transition state” [42].

In Fig. 16, we see an example of the transition state described earlier in this paper. Fixing a periodic Lyapunov orbit about $L_{1}$ (i.e., fixing $I_{2} = 1.0$ ) and fixing a point along the orbit (i.e., fixing $θ_{2} = 0$ ), we generate four sample trajectories, with one in each quadrant of the local saddle manifold, by fixing $I_{1} = 0.125$ and initializing $θ_{1} = 0$ with an imaginary offset for quadrants II–IV. The blue and the orange trajectories, which pass through the transition state clearly in the right-hand-side figure, move from the Earth (Moon) side to the Moon (Earth) side, respectively. On the other hand, the trajectories in the other quadrants (purple and green) do not cross the transition state, and hence remain on the Earth and Moon sides, respectively. Because the CR3BP is a chaotic dynamical system, we expect points in phase space to move in and out of neighborhoods about the equilibria over time, possibly passing through multiple regions. Hence, a trajectory will have different representations while in distinct regions; but using the notion of the transition state, we can recover information about where an unbounded trajectory came from and where it is going.

Figure 17 shows a point near the $L_{1}$ region that moves into the $L_{2}$ region, showing both action-angle element and the corresponding Lissajous orbits. This trajectory is generated by an initial quasi-periodic invariant torus near $L_{1}$ at a sufficiently high energy level, perturbing along the unstable manifold and numerically integrating until the trajectory becomes sufficiently close to $L_{2}$ . Stopping the integration, we transform to action-angle elements about $L_{2}$ to determine the corresponding $L_{2}$ Lissajous invariant torus. Note that if the corresponding saddle integral $I_{1}$ about $L_{2}$ were approximately zero, then the trajectory would have to lie on the stable manifold of some $L_{2}$ Lissajous orbit. This approach was used in Ref. [43]. Moreover, this trajectory lies locally on the $L_{1}$ saddle manifold of the $L_{1}$ Lissajous orbit, which is highlighted on the left, with the black diamond indicating the point along the orbit with the saddle components set to zero. Similarly, this trajectory also lies locally on the $L_{2}$ saddle manifold of the $L_{2}$ Lissajous orbit, highlighted on the right at a later time, with the black diamond indicating the corresponding point as well. Due to this, the local action-angle orbital elements are not necessarily unique across the larger phase space, even though they can be computed locally.

Fig. 17
Local action-angle orbital element corresponding tori near Earth–Moon $L_{1}$ and $L_{2}$ .

VII. Mission Applications

A. Generalizing to the Saddle Manifold

The local action-angle orbital elements allow us to exploit a larger portion of the saddle manifold rather than two of its submanifolds: the stable and unstable manifolds. Using a larger subset of the higher-dimensional saddle manifold offers an additional dimension of options for trajectory design, on top of understanding transiting and nontransiting trajectories through cislunar and translunar regions (more generally, $L_{1}$ and $L_{2}$ in any planet–Moon or Sun–planet CR3BP). We can understand the saddle manifold as the more general “tube” taking a spacecraft from one region to another without using additional propellant. We can see this in Fig. 16. Although the stable manifold will take a spacecraft directly to the highlighted Lyapunov planar periodic orbit, the saddle trajectories come close to the orbit but either transit to/from $L_{2}$ or remain on the $L_{1}$ or $L_{2}$ side (nontransit).

Moreover, the saddle plane allows us to describe the winding conditions of trajectories near the Lagrange points. A smaller value of $I_{1}$ indicates a trajectory closer to the stable or unstable manifold; hence, the desired number of windings about an equilibrium point could be determined by studying an upper or lower bound on $I_{1}$ . Additionally, with the local action-angle orbital element parametrization, the corresponding orbit (in this case, a planar Lyapunov periodic orbit or a Lissajous quasi-periodic in the out-of-plane case) is automatically parametrized by the remaining actions ( $I_{2}$ , $I_{3}$ ).

Understanding the transit/nontransit conditions, as well as the winding condition about a libration point, can ease the design of a mission such as Genesis (which used the unstable manifold of a Sun–Earth $L_{1}$ halo orbit, transited to the Sun–Earth $L_{2}$ side, wound once about $L_{2}$ , and then returned to Earth) with a parking orbit option [29] (see Fig. 18). The local action-angle orbital elements presented can be used between libration points to simplify complex preliminary mission designs. Because we have shown that halo orbits can still be interpreted using local action-angle orbital elements despite being outside the radius of convergence of the complete normal form, for the case of Genesis, in which a halo orbit about $L_{1}$ was used, an action-angle element representation could be used to ease the generation of a transiting trajectory to the $L_{2}$ side before returning to a parking orbit. For design about a particular target periodic orbit, one can modify this procedure to recenter about the target orbit, as discussed in a previous section. In the Genesis example, for the trajectory arc about $L_{1}$ , the local saddle manifold orbital elements can be used to determine the transiting condition: in particular, trajectories in quadrant I transit from the $L_{1}$ side to the $L_{2}$ side. Hence, small positive perturbations of the initial $L_{1}$ orbit onto the unstable manifold, parametrized by $q_{1}$ , will generate several windings before transiting to the $L_{2}$ side. Conversely, around $L_{2}$ , a nontransit condition is used to place the trajectory in quadrant IV of the saddle plane. Note that the quadrants of the saddle plane can be identified by using the normal form variables ( $q_{1}$ , $p_{1}$ ) or by the local action-angle pair ( $I_{1}$ , $φ_{1}$ ) by considering the imaginary offset in the hyperbolic angular variable $φ_{1}$ .

Fig. 18
Genesis mission trajectory and flight plan [29].

B. Targeting

Describing the local dynamics in action-angle coordinates simplifies the equations of motion to constant actions and linearly changing angles. The stable and unstable manifolds are parametrized as $W^{s} = {p_{1} \in R | I_{1} = 0 \Rightarrow q_{1} = 0}$ and $W^{u} = {q_{1} \in R | I_{1} = 0 \Rightarrow p_{1} = 0}$ , respectively, so that $p_{1}$ (respectively, $q_{1}$ ) parametrizes the point along the stable (respectively, unstable) manifold. Simplifying the parametrization of desired trajectories can enable improved targeting performance. For example, stationkeeping using a saddle generalization, or the Floquet mode approach, is simplified by having a local parametrization of trajectories in the vicinity. Because stationkeeping is often a local procedure, having approximately integrable dynamics greatly reduces the complexity of the problem; and the visualization tools described here can be used. Being able to target specific actions, without necessitating the inclusion of the angle variables, would reduce the dimensionality of the problem by half.

In addition to stationkeeping, another example application of local action-angle elements is the transfer problem. To illustrate this further, consider the Artemis-P1 (2007) spacecraft, which performed a transfer from Earth–Moon $L_{2}$ to $L_{1}$ [27]. By using local action-angle elements, finding homoclinic and heteroclinic connections, i.e., free transfers, between orbits in cislunar space reduces to solving a boundary value problem with $I_{1} = 0$ in the region of validity. To target a particular energy level, one can use the energy level sets described in the error analysis section because the energy level set enables searching over a family of quasi-periodic orbits parametrized by the energy level set rather than by computation of individual tori. To target a specific libration point orbit, one can incorporate conditions on the corresponding actions: $I_{2}$ and $I_{3}$ . Although these connections will not be as robust as a numerical differential corrections algorithm, they can be used as an initial guess, which is required by such numerical techniques. Once again, the local action-angle orbital element framework can be used to simplify complex preliminary mission designs such as Artemis-P1.

C. Benefits of the Complete Normal Form

One can perform similar saddle manifold analyses using a seminormal form called the reduction to the center manifold. In this procedure, the saddle integral $I_{1}$ is still defined but is set to zero to parametrize to the center manifold. Although the reduction to the center manifold in general achieves a larger radius of convergence, due to the removal of fewer terms in the Hamiltonian, one of the benefits of using a complete normal form over the reduction to the center manifold is that, in a complete normal form, one is able to immediately identify and parametrize the corresponding orbit or invariant surface in targeting, mission analysis, etc. This parametrization gives the engineer more control due to the ease of generating invariant objects and orbits quickly. Additionally, this parametrization exploits higher-order expansions than previous methods [32,33].

VIII. Data Considerations

The normal forms computed to order 32 (order 16 in the action variables) about the collinear equilibria each contain 968 terms; the direct and inverse coordinate transformation-generating functions are stored, with each using 0.35658 GB. Similarly, the normal forms computed to order 32 (order 16 in the action variables) about the triangular equilibria each contain 968 terms; the direct and inverse coordinate transformation-generating functions are stored each using 0.064836 and 0.16939 GB, respectively. Note that the direct and inverse coordinate transformation generating functions need to be computed only once and stored. Moreover, one can compute transformations at any order that is at or lower than the stored order because the generating functions required for lower-order computation are a subset of the generating functions required for a higher-order computation:

{G_{k}}_{k = 1, \dots, n} \subseteq {G_{k}}_{k = 1, \dots, n, \dots, N}

IX. Conclusions

In this paper, local action-angle orbital elements were constructed, which are local semi-analytical analogs to orbital elements in the circular restricted three-body problem that are defined using action-angle coordinates in the Birkhoff–Gustavson normal form about a bounded special solution. By using local action-angle orbital elements centered around the five equilibria, the dynamics of the region (approximate up to a finite order) can be explicitly solved for, families of periodic and quasi-periodic orbits (in relation to each other) are parametrized, and the elliptic and hyperbolic modes of motion are separated out. By explicitly solving for the dynamics up to a finite order, the nonlinear dynamics become integrable; and the dynamics around the libration points are greatly simplified. Local action-angle orbital elements were shown to provide a practical complement to numerical methods for cataloging periodic and quasi-periodic orbits. Additionally, this high-order parametrization offers a novel exploitation of the entire saddle manifold that is beyond the typical stable and unstable manifolds, which are submanifolds of the saddle. Using action-angle elements to understand the transition state, one can determine if an unbounded trajectory is transit or nontransit; in other words, given a spacecraft state with an unknown time history in the vicinity of a libration point, the imaginary offset in the hyperbolic angle element reduces the set of options from where in the system a spacecraft came. Finally, it was found that the local action-angle elements may have varying definitions for distinct regions of phase space; e.g., transiting trajectories around $L_{1}$ and $L_{2}$ correspond to distinct invariant tori. The implication, however, is that every point in phase space will have at least one local action-angle orbital element representation near some bounded invariant manifold. Moreover, the Genesis and Artemis-P1 trajectories offer further motivation to consider transitioning action-angle elements between $L_{1}$ and $L_{2}$ in both Sun–Earth and Earth–Moon frames.

Appendix: Local Action-Angle Orbital Elements for a Saddle

The purpose of this Appendix is to define the local action-angle orbital elements for a saddle; in particular, these elements are given by either the hyperbolic action-angle coordinates, which are singular along the stable and unstable manifolds, or the standard normal form coordinates corresponding to the saddle.

Consider a two-dimensional saddle, which is given by the position-momenta pair ( $q$ , $p$ ). The equations of motion are given by

\dot{q} = q_{0} e^{λ t}, \dot{p} = p_{0} e^{- λ t}

(A1)

Hence, the solutions of the saddle are organized as follows: $(q, p) = (0, 0)$ is a fixed point; along the $q$ axis, the flow goes toward $\pm \infty$ exponentially in forward time; along the $p$ axis, the flow goes toward $\pm \infty$ exponentially in backward time; and the remaining solutions are hyperbolic trajectories in each quadrant.

In local action-angle orbital elements on a center manifold [i.e., action-angle variables (similar to a harmonic oscillator) near an elliptic fixed point], the action variable fixes a trajectory, and the angle variable indicates a point along the trajectory. For the action-angle elements on a saddle manifold, we follow a similar intuition.

Because the solutions remain in the same quadrant as their initial condition, we will state a definition and check its validity in each quadrant. As discussed previously, and from Eq. (A1), the saddle emits an integral: $I = q p$ . Because we define the action-angle elements for a center manifold by integral and conjugate angles, here, we define the hyperbolic angle coordinate to be the symplectic conjugate coordinate to the saddle integral. Namely, we identify the transformation:

(q, p) \mapsto (I, θ) ≔ (q p, \ln \frac{\sqrt{q}}{\sqrt{p}})

(A2)

Note that $(I, θ) \in R_{\geq 0} \times R$ for ( $q$ , $p$ ) in the first quadrant, i.e., $q$ , $p \geq 0$ . One can easily check that the preceding transformation is symplectic. The inverse transformation is given by

(I, θ) \mapsto (\sqrt{I} e^{θ}, \sqrt{I} e^{- θ}) = (q, p)

(A3)

Furthermore, the equations of motion in the hyperbolic action-angle coordinates are $\dot{I} = 0$ and $\dot{θ} = λ$ .

Now, for the second quadrant (i.e., ${q < 0, p > 0}$ ), the domain of ( $I$ , $θ$ ) becomes complex. Immediately, we see that $I \in R_{\leq 0}$ . To find the domain of $θ$ , we compute

θ = \ln \frac{\sqrt{q}}{\sqrt{p}} = \ln \frac{i \sqrt{| q |}}{\sqrt{| p |}} = \ln (\frac{\sqrt{| q |}}{\sqrt{| p |}}) + i π / 2

(A4)

Hence, $(I, θ) \in R_{\leq 0} \times R + i π / 2$ in the second quadrant. We check the composition of direct and inverse transformations to see that we obtain $(q, p) \mapsto (I, θ) \mapsto (q, p)$ :

q = \sqrt{I} e^{θ} = i \sqrt{| q | | p |} e^{\ln (\sqrt{| q |} / \sqrt{| p |}) + i π / 2} = i^{2} | q | = - | q |

(A5)

p = \sqrt{I} e^{- θ} = i \sqrt{| q | | p |} e^{- \ln (\sqrt{| q |} / \sqrt{| p |}) - i π / 2} = - i^{2} | p | = | p |

(A6)

as desired because

q = - | q |

and

p = | p |

in the second quadrant.

For the third quadrant (i.e., { $q$ , $p < 0$ }), we perform the composite direct and inverse transformation to see the imaginary offset for the hyperbolic angle coordinate:

q = I e^{θ} = I e^{R e (θ) + i Im (θ)} = \sqrt{q p} \sqrt{\frac{| q |}{| p |}} e^{i Im (θ)} = | q | e^{i Im (θ)}

(A7)

Because $q = - | q |$ in the third quadrant, we see that $e^{i Im (θ)} = - 1$ , and so $Im (θ) = π$ . In this way, we obtain $(I, θ) \in R_{\geq 0} \times R + i π$ .

In the fourth quadrant, a similar computation to quadrant 2 shows $(I, θ) \in R_{\leq 0} \times R + i 3 π / 2$ .

1. Singularity Along $q$ and $p$ Axes

Observe along the $q$ and $p$ axes that $I = 0$ and $θ = \mp \infty$ , respectively, where the sign indicates the $q$ or $p$ axis (i.e., the stable or unstable manifold) and where $θ$ has some imaginary offset of either $i 0$ , $i π / 2$ , $i π$ or $i 3 π / 2$ , depending on the quadrant. The four imaginary offsets will be paired corresponding to the positive or negative stable/unstable manifold. Accordingly, the angle coordinate needs infinite time to move between points along the same axis; hence, this is a singularity of the hyperbolic action-angle coordinates on the saddle. Because the singularity occurs on a set with a measure of zero in the phase space, this is reminiscent of the singularity in the Keplerian orbital elements.

2. Application for Nonsingular Local Action-Angle Orbital Elements

On a saddle manifold, we define local action-angle orbital elements to be the action-angle pair ( $I$ , $θ$ ), where the imaginary offset of $θ$ indicates the quadrant of motion. Note that we could also check the signs of $q$ and $p$ to determine the quadrant and use the formulas accordingly. We then use an abuse of notation, where we take the absolute value of the integral $| I |$ and the real part of $θ$ [ $R e (θ)$ _] as the coordinates ( $I$ , $θ$ ). Now, if $| I | = 0$ or, equivalently, if $sgn (q) = 0$ or $sgn (p) = 0$ , then the coordinates are singular. In this case, we retain the ( $q$ , $p$ ) normal form coordinates as the action-angle elements. These act as action-angle coordinates, similar to the ( $I$ , $θ$ ), wherein $I$ fixes a trajectory and $θ$ fixes a point along the trajectory. Along the stable or unstable manifold, similarly, if $q = 0$ (without loss of generality), then the stable manifold is the fixed trajectory; and the value of $p$ fixes a point along the stable manifold. Note that $sgn (p)$ [ $sgn (q)$ , respectively] determines the positive or negative stable manifold. By using hyperbolic action-angle coordinates with the normal form coordinates on the singularity, we define a set of nonsingular local action-angle orbital elements on a saddle.

Acknowledgments

This work was supported by a National Defense Science & Engineering Graduate Fellowship. Daniel J. Scheeres acknowledges support from U.S. Air Force Office of Scientific Research grant FA9550-21-1-0332. This work was presented as paper American Astronautical Society 22-756 at the 2022 AAS/AIAA Astrodynamics Specialist Conference in August 2022.

References

[1] Roy A. E., Orbital Motion, CRC Press, Boca Raton, FL, 2020, Chap. 4. https://doi.org/10.1201/9780367806620 Google Scholar
[2] Meyer K. R. and Hall G. R., Introduction to Hamiltonian Dynamical Systems and the N-Body Problem, Springer-Verlag, New York, 1992, Chap. 3. https://doi.org/10.1007/978-1-4757-4073-8_1 Google Scholar
[3] Scheeres D. J., Orbital Motion in Strongly Perturbed Environments, Springer, Berlin, 2012, Chap. 4. https://doi.org/10.1007/978-3-642-03256-1 Crossref Google Scholar
[4] Delaunay C., Théorie du Mouvement de la Lune, Vol. 28, Mallet-Bachelier, 1867, pp. 1–15. Google Scholar
[5] Gurfil P. and Kasdin N. J., “Canonical Modelling of Coorbital Motion in Hill’s Problem Using Epicyclic Orbital Elements,” Astronomy and Astrophysics, Vol. 409, No. 3, 2003, pp. 1135–1140. https://doi.org/10.1051/0004-6361:20031162 Crossref Google Scholar
[6] Danby J., Fundamentals of Celestial Mechanics, Willman-Bell, Richmond, 1992, Chap. 8. https://doi.org/10.1007/bf00051051 Google Scholar
[7] Arnold V. I., Mathematical Methods of Classical Mechanics, Graduate Texts in Mathematics, Vol. 60, Springer, New York, 1978, Chap. 2. https://doi.org/10.1007/978-1-4757-1693-1 Google Scholar
[8] Birkhoff G. D., Dynamical Systems, Vol. 9, American Mathematical Soc., Providence, RI, 1927. https://doi.org/10.1090/coll/009 Google Scholar
[9] Arnold V. I., Kozlov V. V. and Neishtadt A. I., Mathematical Aspects of Classical and Celestial Mechanics, 3rd ed., Encyclopaedia of Mathematical Sciences Series, Springer, Berlin, 2006, Appendix 7. https://doi.org/10.1007/978-3-540-48926-9 Crossref Google Scholar
[10] Moser J., “New Aspects in the Theory of Stability of Hamiltonian Systems,” Communications on Pure and Applied Mathematics, Vol. 11, No. 1, 1958, pp. 81–114. https://doi.org/10.1002/cpa.3160110105 Crossref Google Scholar
[11] Jorba À. and Masdemont J., “Dynamics in the Center Manifold of the Collinear Points of the Restricted Three Body Problem,” Physica D: Nonlinear Phenomena, Vol. 132, Nos. 1–2, 1999, pp. 189–213. https://doi.org/10.1016/s0167-2789(99)00042-1 Crossref Google Scholar
[12] Poincaré H., New Methods of Celestial Mechanics: Volume 1 Periodic Solutions. Non-Existence of Uniform Integrals. Asymptotic Solutions, Dover, New York, 1957. https://doi.org/10.5860/choice.30-4974 Google Scholar
[13] Gustavson F., “On Constructing Formal Integrals of a Hamiltonian System Near an Equilibrium Point,” Astronomical Journal, Vol. 71, Oct. 1966, Paper 670. https://doi.org/10.1086/110172 Google Scholar
[14] Von Zeipel H., “Recherches sur le mouvement des petites planetes,” Arkiv för Matematik, Astronomi och Fysik, Vol. 11, Nos. 1–2, 1916, pp. 1–58. Google Scholar
[15] Giorgilli A. and Galgani L., “Formal Integrals for an Autonomous Hamiltonian System Near an Equilibrium Point,” Celestial Mechanics, Vol. 17, No. 3, 1978, pp. 267–280. https://doi.org/10.1007/bf01232832 Google Scholar
[16] Deprit A., “Canonical Transformations Depending on a Small Parameter,” Celestial Mechanics, Vol. 1, No. 1, 1969, pp. 12–30. https://doi.org/10.1007/bf01230629 Crossref Google Scholar
[17] Jorba À., “A Methodology for the Numerical Computation of Normal Forms, Centre Manifolds and First Integrals of Hamiltonian Systems,” Experimental Mathematics, Vol. 8, No. 2, 1999, pp. 155–195. https://doi.org/10.1080/10586458.1999.10504397 Crossref Google Scholar
[18] Jorba A. and Villanueva J., “Numerical Computation of Normal Forms Around Some Periodic Orbits of the Restricted Three-Body Problem,” Physica D: Nonlinear Phenomena, Vol. 114, Nos. 3–4, 1998, pp. 197–229. https://doi.org/10.1016/s0167-2789(97)00194-2 Google Scholar
[19] Subramanian S. C., Waswa P. M. and Redkar S., “Computation of Lyapunov–Perron Transformation for Linear Quasi-Periodic Systems,” Journal of Vibration and Control, Vol. 28, Nos. 11–12, 2022, pp. 1402–1417. https://doi.org/10.1177/1077546321993568 Google Scholar
[20] Duarte Ferreira G., “On the Dynamics Around the Collinear Points in the Sun-Jupiter System,” Ph.D. Thesis, Univ. of Barcelona, Barcelona, Spain, 2020. Google Scholar
[21] Whitley R. and Martinez R., “Options for Staging Orbits in Cislunar Space,” 2016 IEEE Aerospace Conference, IEEE, New York, 2016, pp. 1–9. https://doi.org/10.1109/aero.2016.7500635 Google Scholar
[22] Crusan J., Bleacher J., Caram J., Craig D., Goodliff K., Herrmann N., Mahoney E. and Smith M., “NASA’s Gateway: An Update on Progress and Plans for Extending Human Presence to Cislunar Space,” 2019 IEEE Aerospace Conference, IEEE, New York, 2019, pp. 1–19. https://doi.org/10.1109/aero.2019.8741561 Google Scholar
[23] Holzinger M., Chow C. and Garretson P., “A Primer on Cislunar Space,” U.S. Air Force Research Lab. Space Vehicles Directorate 2021-1271, 2021, pp. 1–23, https://www.afrl.af.mil/Portals/90/Documents/RV/A%20Primer%20on%20Cislunar%20Space_Dist%20A_PA2021-1271.pdf?ver=vs6e0sE4PuJ51QC-15DEfg%3D%3D. Google Scholar
[24] Davis D. C., Phillips S. M., Howell K. C., Vutukuri S. and McCarthy B. P., “Stationkeeping and Transfer Trajectory Design for Spacecraft in Cislunar Space,” AAS/AIAA Astrodynamics Specialist Conference, Springer Nature, London, 2017, pp. 1–20. Google Scholar
[25] Strange N., Landau D., McElrath T., Lantoine G. and Lam T., “Overview of Mission Design for NASA Asteroid Redirect Robotic Mission Concept,” 33rd International Electric Propulsion Conference (IEPC2013), Electric Rocket Propulsion Soc., Brookpark, OH, 2013, Paper 4436. Google Scholar
[26] Parrish N. L. O., “Low Thrust Trajectory Optimization in Cislunar and Translunar Space,” Ph.D. Thesis, Univ. of Colorado Boulder, Boulder, CO, 2018. Google Scholar
[27] Folta D. C., Woodard M., Howell K., Patterson C. and Schlei W., “Applications of Multi-Body Dynamical Environments: The ARTEMIS Transfer Trajectory Design,” Acta Astronautica, Vol. 73, April 2012, pp. 237–249. https://doi.org/10.1016/j.actaastro.2011.11.007 Crossref Google Scholar
[28] Folta D. C., Pavlak T. A., Haapala A. F., Howell K. C. and Woodard M. A., “Earth–Moon Libration Point Orbit Stationkeeping: Theory, Modeling, and Operations,” Acta Astronautica, Vol. 94, No. 1, 2014, pp. 421–433. https://doi.org/10.1016/j.actaastro.2013.01.022 Crossref Google Scholar
[29] Lo M. W., Williams B. G., Bollman W. E., Han D., Hahn Y., Bell J. L., Hirst E. A., Corwin R. A., Hong P. E., Howell K. C. and Barden B., “Genesis Mission Design,” Journal of the Astronautical Sciences, Vol. 49, No. 1, 2001, pp. 169–184. https://doi.org/10.1007/bf03546342 Crossref Google Scholar
[30] Folta D. C., Bosanac N., Guzzetti D. and Howell K. C., “An Earth–Moon System Trajectory Design Reference Catalog,” Acta Astronautica, Vol. 110, May 2015, pp. 341–353. https://doi.org/10.1016/j.actaastro.2014.07.037 Crossref Google Scholar
[31] Guzzetti D., Bosanac N., Haapala A., Howell K. C. and Folta D. C., “Rapid Trajectory Design in the Earth–Moon Ephemeris System via an Interactive Catalog of Periodic and Quasi-Periodic Orbits,” Acta Astronautica, Vol. 126, Sept. 2016, pp. 439–455. https://doi.org/10.1016/j.actaastro.2016.06.029 Crossref Google Scholar
[32] Burnett E. R. and Schaub H., “Spacecraft Relative Motion Dynamics and Control Using Fundamental Modal Solution Constants,” Journal of Guidance, Control, and Dynamics, Vol. 45, No. 10, 2022, pp. 1786–1799. https://doi.org/10.2514/1.g006603 Link Google Scholar
[33] Elliott I. and Bosanac N., “Describing Relative Motion Near Periodic Orbits via Local Toroidal Coordinates,” Celestial Mechanics and Dynamical Astronomy, Vol. 134, No. 2, 2022, Paper 19. https://doi.org/10.1007/s10569-022-10074-8 Google Scholar
[34] Mireles James J. D. and Murray M., “Chebyshev–Taylor Parameterization of Stable/Unstable Manifolds for Periodic Orbits: Implementation and Applications,” International Journal of Bifurcation and Chaos, Vol. 27, No. 14, 2017, Paper 1730050. https://doi.org/10.1142/s0218127417300506 Google Scholar
[35] de la Llave R., Marco J. M. and Moriyón R., “Canonical Perturbation Theory of Anosov Systems and Regularity Results for the Livsic Cohomology Equation,” Annals of Mathematics, Vol. 123, No. 3, 1986, pp. 537–611. https://doi.org/10.2307/1971334 Google Scholar
[36] Szebehely V., “Modifications of the Restricted Problem,” Theory of Orbit, Elsevier, New York, 1967, pp. 556–652. https://doi.org/10.1016/B978-0-12-395732-0.50016-7 Google Scholar
[37] Jorba À., Jorba-Cuscó M. and Rosales J. J., “The Vicinity of the Earth–Moon L1 Point in the Bicircular Problem,” Celestial Mechanics and Dynamical Astronomy, Vol. 132, No. 2, 2020, p. 11. https://doi.org/10.1007/s10569-019-9940-2 Google Scholar
[38] Duarte G. and Jorba À., “Using Normal Forms to Study Oterma’s Transition in the Planar RTBP,” Discrete and Continuous Dynamical Systems-B, Vol. 28, No. 1, 2022, pp. 230–244. https://doi.org/10.3934/dcdsb.2022073 Google Scholar
[39] Giorgilli A., Notes on Hamiltonian Dynamical Systems, Vol. 102, Cambridge Univ. Press, Cambridge, England, U.K., 2022, Chap. 5. https://doi.org/10.1017/9781009151122 Google Scholar
[40] Paez R. I. and Guzzo M., “On the Semi-Analytical Construction of Halo Orbits and Halo Tubes in the Elliptic Restricted Three-Body Problem,” Physica D: Nonlinear Phenomena, Vol. 43, Nov. 2022, Paper 133402. https://doi.org/10.1016/j.physd.2022.133402 Google Scholar
[41] Baresi N., Olikara Z. P. and Scheeres D. J., “Fully Numerical Methods for Continuing Families of Quasi-Periodic Invariant Tori in Astrodynamics,” Journal of the Astronautical Sciences, Vol. 65, No. 2, 2018, pp. 157–182. https://doi.org/10.1007/s40295-017-0124-6 Crossref Google Scholar
[42] Uzer T., Jaffé C., Palacián J., Yanguas P. and Wiggins S., “The Geometry of Reaction Dynamics,” Nonlinearity, Vol. 15, No. 4, 2002, Paper 957. https://doi.org/10.1088/0951-7715/15/4/301 Google Scholar
[43] Waalkens H., Burbanks A. and Wiggins S., “A Computational Procedure to Detect a New Type of High-Dimensional Chaotic Saddle and Its Application to the 3D Hill’s Problem,” Journal of Physics A: Mathematical and General, Vol. 37, No. 24, 2004, Paper L257. https://doi.org/10.1088/0305-4470/37/24/l04 Google Scholar

Tables

**Table 1 Differences between a normal form prediction and a numerical integration about $L_{1}$**
$λ_{0}$	$‖ v_{1} - v_{0}^{1} ‖_{2}$
0.00001	2.985839760499337e-08
0.00002	8.445820541952098e-06
0.00004	1.193880187497116e-05
0.00008	1.687316131245734e-05
0.00016	2.384034784180320e-05
0.00032	3.367097716230715e-05

Journal of Guidance, Control, and Dynamics, volume 46 issue 12 cover

Volume 46, Number 12December 2023

Metrics

Crossmark

Information

Copyright © 2023 by the American Institute of Aeronautics and Astronautics, Inc. All rights reserved. All requests for copying and permission to reprint should be submitted to CCC at www.copyright.com; employ the eISSN 1533-3884 to initiate your request. See also AIAA Rights and Permissions www.aiaa.org/randp.

Topics

Keywords

Acknowledgments

Digital

Received16 December 2022
Accepted17 June 2023
Published online8 August 2023

Local Orbital Elements for the Circular Restricted Three-Body Problem

Abstract

I. Introduction

II. Dynamical Model: The Circular Restricted Three-Body Problem

A. Recentering About Equilibria

B. Series Expansion of Hamiltonian

III. Local Action-Angle Orbital Elements

A. Preliminaries

B. Lie Series Method

C. Local Action-Angle Orbital Elements

D. High-Level Overview

IV. Error Analysis

V. Application to Triangular Equilibria

VI. Application to Collinear Equilibria

VII. Mission Applications

A. Generalizing to the Saddle Manifold

B. Targeting

C. Benefits of the Complete Normal Form

VIII. Data Considerations

IX. Conclusions

Appendix: Local Action-Angle Orbital Elements for a Saddle

1. Singularity Along q and p Axes

2. Application for Nonsingular Local Action-Angle Orbital Elements

Acknowledgments

References

Tables

What's Popular

Metrics

Information

1. Singularity Along $q$ and $p$ Axes