Dynamics Around the Earth-Moon Triangular Points in the Hill Restricted 4-Body Problem

This paper investigates the motion of a small particle moving near the triangular points of the Earth-Moon system. The dynamics are modeled in the Hill restricted 4-body problem (HR4BP), which includes the effect of the Earth and Moon as in the circular restricted 3-body problem (CR3BP), as well as the direct and indirect effect of the Sun as a periodic time-dependent perturbation of the CR3BP. Due to the periodic perturbation, the triangular points of the CR3BP are no longer equilibrium solutions; rather, the triangular points are replaced by periodic orbits with the same period as the perturbation. Additionally, there is a 2:1 resonant periodic orbit that persists from the CR3BP into the HR4BP. In this work, we investigate the dynamics around these invariant objects by computing families of 2-dimensional invariant tori and their linear normal behavior. We identify bifurcations and relationships between families. Mechanisms for transport between Earth, L4, and Moon are discussed. Comparisons are made between the results presented here and in the bicircular problem (BCP).


Introduction
The study of dynamics in the Earth-Moon (EM) environment has become the focus of much research due to the announcement of the lunar gateway placed in a southern L 2 near rectilinear halo orbit (NRHO) [1], as well as nearby orbits used for Artemis program operations [2].As the traffic of spacecraft increases for operations near the Moon, there is an increased need for relay communications between Earth and the vicinity of the Moon.TYCHO is a proposed mission to provide such communication relay services via satellites placed near the Earth-Moon L 4 libration point-one of the two so-called triangular points [3].To model the dynamics in cislunar space, the circular restricted 3-body problem (CR3BP) is used as a first approximation [4,5].In the Earth-Moon CR3BP, the triangular points L 4,5 are elliptic equilibria-they are linearly stable.In fact, it has been shown that there is a region of effective stability around L 4 in the CR3BP [6].This stability, along with the fact that L 4 is equidistant to the Earth and Moon, presents an ideal region of cislunar space to put communications satellites.
The presence of the Trojan asteroids near the Sun-Jupiter triangular points has raised questions about the existence of similar objects in the Earth-Moon system [7].After all, the triangular points are elliptic in both the Sun-Jupiter and Earth-Moon CR3BP, and Earth-Moon Trojans could pose a risk of collision with communications satellites.While it is known that there are no Earth-Moon Trojans, the existence of Kordylewski dust clouds around EM L 4 has recently been confirmed [8,9].The key difference between the Sun-Jupiter and Earth-Moon systems is that there are no other major perturbing forces in the Sun-Jupiter system, as the Sun and Jupiter are the two most massive celestial bodies in the solar system, while the Earth-Moon CR3BP neglects the significant perturbing force of the Sun's gravity.It has been shown that accounting for the gravitational effect of the Sun in the Earth-Moon system qualitatively changes the behavior of L 4 -its stability changes from elliptic to partially hyperbolic [10,11].Accounting for the Sun's gravity helps to describe the non-existence of Earth-Moon Trojan asteroids, and recent work has incorporated additional nongravitational perturbations caused by the Sun to help describe the existence of the Kordylewski clouds from a dynamical astronomy perspective [7,12].However, what are the implications of this qualitative change in dynamical behavior around EM L 4 for communications satellites supporting cislunar operations?The present work seeks to address this question.
The simplest model introducing the periodic forcing of the Sun in the Earth-Moon system is the bicircular restricted 4-body problem (BCP).This model is formed by assuming that the Earth-Moon barycenter travels around the Sun in a circular orbit at a constant rate, while the Earth and Moon move around their barycenter in a circular orbit.The BCP accounts for the Sun's gravitational effect on the particle (the direct effect) but neglects the gravitational effect on the relative motion of the Earth and Moon (the indirect effect).Due to this, the Sun, Earth, and Moon do not form a solution to the 3-body problem, hence we say the model is incoherent.Despite its incoherence, the BCP captures the qualitative behavior of the periodic solutions around L 4 similarly to its coherent counterpart, the quasi-bicircular restricted 4-body problem (QBCP), as shown in Figure 25, borrowed from Andreu [13].There has been some investigation into the dynamics around the triangular points in the BCP, namely the continuation of the L 4 point into its periodic counterparts [11], as well as the computation of quasi-periodic invariant 2-tori in the vertical direction [14].Moreover, Jorba investigated the existence of stable motions near the triangular points of the Earth-Moon system in a high-fidelity ephemeris model, relating to work done in the BCP [15].By a similar group of researchers, there has been more investigation modeling the dynamics affecting the Kordylewski clouds [7].Along the lines of communications satellites being placed near L 4 , some authors have used the BCP as their dynamical model for transfer design [16,17].
In this work we use the Hill restricted 4-body problem (HR4BP) to model the motion of a particle in the Sun-Earth-Moon (SEM) system.The HR4BP is a coherent periodic Hamiltonian system that is both a generalization of the Earth-Moon CR3BP and Sun-Earth Hill's problem [10,18].The Sun-Earth-Moon HR4BP is coherent and models both the direct and indirect effect of the Sun, using a Hill's variational orbit as the relative motion of the Earth and Moon to capture the indirect effect.The HR4BP has been studied minimally in comparison to the BCP and QBCP [10,[18][19][20][21], and the purpose of employing this model in the present work is twofold.On one hand, only the original Scheeres paper which introduced the HR4BP has an investigation of the dynamics around EM L 4 in this model, and even that analysis is austere.So, we would like to further our understanding of the dynamics in this region of the phase space of the system.On the other hand, we would like to compare our results where applicable with the work previously done in the BCP and extend computations to include more than what has been done in the BCP to date.Additionally, from an astrodynamics perspective, transport between Earth and the Moon is of particular interest, so finding any avenues that exist passing through the triangular points may have value for transferring a satellite between cislunar regions of interest.Furthermore, practical stability regions near the triangular points have positive implications for the observability of spacecraft in cislunar space.
The paper is organized as follows.In Section 2 we review the dynamical model used in our investigation, including the continuation of equilibria and resonant periodic orbits near EM L 4 .In Section 3 we give high-level overviews of the several techniques employed to study dynamics around the triangular points.Section 4 makes up the majority of the text, as we detail our findings in the form of families of invariant tori, a Poincaré section, and transport via hyperbolic invariant manifolds.Finally, we draw comparisons where applicable between the existing work in the BCP and our examination of the HR4BP in Section 5.

Hill Restricted 4-Body Problem
The HR4BP is a coherent non-autonomous π-periodic generalization of the CR3BP and the Hill restricted 3-body problem.The derivation of the HR4BP is fully detailed in [10].There are three primary assumptions made in the HR4BP: M 0 ≫ M 1 , M 2 , where M i is the mass of each primary body; the mass of the particle is negligible, i.e., M 3 ≪ 1; and the distance between M 1 and M 2 is much less than their distance to M 0 .In the SEM system, these assumptions are valid, where M 0 is the Sun, M 1 the Earth, and M 2 the Moon.Under these assumptions, the Sun and Earth-Moon barycenter move in a circular orbit and the relative motion of the Earth and Moon follows a specific solution to the Hill 3-body problem, a member of a 1-parameter family of planar periodic orbits known as Hill's variational orbits.
These periodic orbits depend on a parameter, m, which is the synodic period of the Earth-Moon orbit in years.Note for the SEM system, we take m = 0.0808.In this work, we take the reference frame to be the average co-rotating frame of the smaller two primaries, rotating at a constant rate: 1 + 1/m, where 1 is the normalized rotation rate of the Sun and Earth-Moon barycenter frame, and 1/m is the average rotation rate of the Earth-Moon orbit.In addition to the parameter m, the HR4BP depends on a second parameter, µ, the mass ratio of the smaller two primaries.Additionally, we take the following SEM unit normalizations.We normalize the length in the Earth-Moon frame by the mean distance between the Earth and the Moon (384,400 km).We normalize time τ such that 2π time units correspond to a synodic month (29.53 days), similar to [20].
To write the equations of motion, we define the conjugate momenta as p x = ẋ − (1 + m)y, p y = ẏ + (1 + m)x, and p z = ż.The motion of the infinitesimal particle in the Earth-Moon rotating frame is then described as Hamiltonian dynamical system given by: where V (x, y, z, τ ; µ, m) is given by: with r 1−µ and r µ given by: where ξ and η are the normalized Hill variation orbit [22]: where the coefficients a n (m) can be computed following [22]; moreover, we use the coefficients given in [20].Note that the equations of motion are π-periodic in τ , which corresponds to synodic half-month periodicity.Recall that T = 2π ω S , so, as the model is π-periodic, we say that ω S = 2.As a 0 (m) = m 2/3 1 − 2 3 m + O(m 2 ) , we note that if m → 0, the Hamiltonian of the restricted three-body problem is recovered.
We will use both the Hamiltonian and Newtonian formulations in this work, so we present the equations of motion in a Newtonian way: where V is as defined in the Hamiltonian formulation above.Finally, we mention several symmetries present in this model.Using the Hill approximation, the perturbing effect of the Sun is symmetric across the Earth-Moon barycenter, so the equations of motion are invariant under the transformation τ → τ + π.Additionally, the Earth-Moon orbit is π-periodic in this frame, so the system is π-periodic.The equations of motion are unchanged under the transformations: S z : (x, y, z, ẋ, ẏ, ż, τ ) → (x, y, −z, ẋ, ẏ, − ż, τ ).
Thus, it suffices to consider the dynamics around L 4 , as the dynamics around L 5 are obtained by applying the S y symmetry.Observe that this similar symmetry similarly holds in the CR3BP.

Continuation of Equilibria and Resonant Periodic Orbits
The lowest dimensional invariant manifolds in the CR3BP are equilibrium points.We introduce the π-periodic forcing of the Sun by incrementally increasing m from m = 0 to m = 0.0808, i.e., from the EM CR3BP to the SEM HR4BP.Once the periodic forcing is added, the equilibria become π-periodic orbits in general.Of course, there may be bifurcations of this periodic orbit in the continuation process, as seen around L 2 in the HR4BP [19,21].We initialize the continuation of the L 4 point into its π-periodic orbit replacement taking x 0 = x L4 , i.e., the initial state at τ = 0 is taken as the CR3BP L 4 point.Using a standard single shooting algorithm provides the differential corrected periodic orbit [23].Periodic orbits replacing the CR3BP equilibria are not the only periodic solutions of the periodically-forced system, as CR3BP periodic orbits of a resonant period may persist.Whereas the initial point at τ = 0 is taken as the equilibrium point, we must carefully select the initial point for the continuation of a CR3BP periodic orbit into the HR4BP.In the SEM HR4BP, there is a 2:1 resonant periodic orbit which can be continued from the EM CR3BP.The 2:1 resonant orbit comes from the 2π-periodic orbit of the L 4 short-period family in the CR3BP.We utilize sub-harmonic Melnikov theory to predict the points that persist [24].We utilize the energy principle-the change in energy from the start and end of a periodic orbit must be zero-following Cenedese and Haller [25], to construct a Melnikov function necessary for the subharmonic analysis at hand.
To introduce the Melnikov used in this analysis, we write the HR4BP in the Newtonian formulation as r = f (X) + εg(X, τ ), where f (X) is the unperturbed system, g(X, τ ) is the π-periodic perturbation, and ε > 0 is small.Expanding the HR4BP equations of motion about m = 0, we have εg = mg 1 + m 2 g 2 + m 3 g 3 + • • • .Brown et al. articulated that g 1 ≡ 0, so we use g 2 as the next lowest order perturbation [26].Repeating the computation of Brown et al., and using α := τ 0 + τ , we have: where R 1−µ,C and R µ,C are the relative positions of the particle with respect to the Earth and the Moon, respectively, in the CR3BP.Originally derived in Brown et al. [26], the Melnikov function for the 2:1 resonant periodic orbit in the HR4BP then takes the form: where X(s + τ ; ε = 0) is the state along the orbit (where s parameterizes the initial condition at τ 0 ) in the unperturbed problem.The simple zeros of M represent the initial points along the orbit where no work is done by non-conservative forces over the orbit-a criterion for periodicity in the non-autonomous system.The top plot in Figure 1 shows the Melnikov function evaluated at each point along the 2:1 resonant orbit in the unperturbed problem.As there are 4 zeros, we continue the periodic orbit into the HR4BP from these points-observe that they are evenly spaced along the orbit and that opposite points correspond to the same periodic orbit with a phase shift.The bottom plot in Figure 1 shows the continuation in m from the EM CR3BP (m = 0) into the SEM HR4BP (m = 0.0808).Notice that there is a period-doubling bifurcation between Branches A and B which destroys Branch B. Hence, only Branches A and C persist in the SEM system.Figure 2 shows the surviving periodic orbits in the configuration space.We call Branch A the "dynamical equivalent of L 4 " or DE L 4 for short; we call Branch C the 2:1 resonant orbit.
The periodic forcing of the Sun qualitatively changes the linear normal behavior of L 4 .So, we present in Table 1 the stability data of the periodic orbits in the SEM HR4BP, i.e., the eigenvalues of the monodromy matrices, as well as the stability of the L 4 CR3BP equilibrium point and 2πperiodic orbit.The π-stroboscopic map is symplectic so the monodromy matrices are symplectic, hence the eigenvalues come in reciprocal pairs.Observe that the L 4 equilibrium point in the EM CR3BP is elliptic.In the SEM HR4BP, the elliptic L 4 point is replaced by a partially hyperbolic DE L 4 π-periodic orbit, i.e., having stability type center × saddle × center.This means there is a 1-parameter family of 2-dimensional invariant tori in the planar and vertical directions, and by coupling oscillations, there is a 2-parameter family of 3-dimensional tori.Conversely, the 2π-periodic orbit in the EM CR3BP is elliptic and continues into an elliptic periodic orbit in the SEM HR4BP, though picks up a center mode because the unity eigenvalues transition into an additional center mode due to the periodic perturbation.This orbit-Branch C-is elliptic over the m considered.Branch B, stemming from the same CR3BP periodic orbit, is partially hyperbolic throughout m for which it persists.Branch B is killed in the bifurcation with Branch A, whence Branch A changes stability from elliptic to partially hyperbolic, inheriting the (un)stable manifold structure through the bifurcation.Figure 3 shows the π-stroboscopic map of the (un)stable manifolds of Branches A and B for several values of m, namely before and after the bifurcation.In this figure, we see the structure of the (un)stable manifolds persists from Branch B to Branch A through the bifurcation.Note that these are not homoclinic points of the flow map, i.e., the (un)stable manifolds do not intersect in phase space.In the SEM HR4BP, Scheeres showed that the stable manifold of DE L 4 is confined to the interior of the structure, and the unstable manifold remains on the exterior.

Methods
In this section, we review the myriad of methods employed to study the dynamics near the Earth-Moon triangular points in the HR4BP.Firstly, we would like to characterize the bounded motions around DE L 4 .Since DE L 4 is partially hyperbolic, a common strategy is to apply center manifold reduction wherein the bounded and hyperbolic motions are systematically decoupled up to some order of a series expansion around DE L 4 .The primary limitation of the center manifold reduction is the radius of convergence which relates to the distance the expansion accurately represents the dynamics around DE L 4 .Oftentimes center manifold reduction can be used in conjunction with the computation of quasi-periodic invariant tori to provide a more complete picture of the dynamics in a particular region; however, if the center manifold reduction technique is severely limited, this can give some insight into the dynamics qualitatively because there is a plausible explanation for this.Hence, we also apply the flow map method [5,27] to compute invariant 2-tori around DE L 4 and the 2:1 resonant orbit.We are also interested in the unbounded motions around DE L 4 , so we compute the normal bundles of invariant tori.Using the stable and unstable bundles, we compute the global stable and unstable manifolds of particular 2-tori.

Center Manifold Reduction
The normal form computational techniques used here are based on [28,29] and modified to the HR4BP, as in [18].We will give an overview of the fundamental aspects of the techniques; see the referenced texts for details on implementation.In this respect, the primary objective of this normal form technique is to re-center the Hamiltonian function about the dynamical equivalents of L1,2, expand as an infinite series, and, through successive time-dependent canonical transformations, change coordinates such that the system is autonomous and the local hyperbolic and elliptic modes of motion are decoupled.See the appendices of [18] for the derivations used for re-centering and expanding the Hamiltonian around EM L 1,2 -the procedure for EM L 4 is identical.
We re-center and expand the Hamiltonian in Equation 1 about Earth-Moon L 4 , denoted by H EM .See [18] for details.This cancels terms of degree one in the Hamiltonian.In addition to recentering, it is useful to autonomize the Hamiltonian by introducing a variable that is a symplectic conjugate to the periodic variable τ , called I S .Then, we define the autonomized Hamiltonian as: Note that the autonomization procedure above does not remove the dependence on τ .We would like to remove the τ -dependence up to a certain order, and this is reserved for a later step.
Once the Hamiltonian function has been autonomized by introducing I S conjugate to τ and re-centered about the dynamical equivalent of interest to cancel terms of degree one, we utilize the Symplectic Floquet Theorem [30] to put the quadratic terms into diagonal normal form with constant coefficients.Note that for degrees one and two, we need not expand the Hamiltonian as an infinite series of homogeneous polynomials [31].Applying the Floquet transformation to the Hamiltonian in Equation 9, as well as the complexification given by for j = 2, 3, we obtain the desired form: We now follow a similar discussion as [28] to put the third-order terms into the desired form.The procedure is then applied similarly for subsequent orders.As mentioned above, we start with the Hamiltonian in Equation 11, which is of the form: , where H 2 is in complex normal form and H n = H n (q, p, τ ), n ≥ 3, is a homogeneous polynomial of degree n whose coefficients are complex-valued periodic functions of τ .
The normalization process begins at degree 3, using a generating function G 3 which is also a homogeneous polynomial of degree 3 in (q, p) with coefficients depending periodically on τ .Then, we continue at degree 4 with G 4 , and so on.By choosing a generating function that defines a non-autonomous Hamiltonian vector field, and flowing along for unit time, we can define a timedependent canonical transformation that removes the time dependence of the Hamiltonian in the new coordinates.
Starting with the following Hamiltonian: where , is a multi-index.Note that we use Fourier series to approximate the time-periodic coefficients a k n (τ ).To remove certain terms from H 3 (q, p, τ ) (similarly for H n (q, p, τ ) for each n), we will make a change of variables generated by G 3 : As G 3 generates a Hamiltonian vector field, the canonical transformation generated by G 3 is the time one flow along this vector field.We write explicitly the terms of degree 3 of the transformed Hamiltonian H3 : where we define the vector of coefficients of H 2 as ω = (ω 1 , iω 2 , iω 3 ).Note that we consider here the stability case of saddle × center × center, but this procedure is easily generalized to other stability cases.Continuing, note that H3 and G 3 are unknowns in Equation 14, and we solve for G 3 presupposing H3 to have a desired form.Imposing the condition that the form of H3 = |k|=3 h k 3 (τ )q k 1 p k 2 , and grouping all terms with the same k, we obtain the following set of linear differential equations with g k 3 as unknowns to be selected later: We can solve the above differential equations explicitly as: where is the "resonance module," i.e., the set of indices of resonant terms.In the Sun-Earth-Moon HR4BP, J k contains only the zero vector for each k.Note that if j ∈ J k , we have to impose a k 3,j = h k 3,j and are unable to remove the time dependence.But, for j / ∈ J k , we can choose h k 3,j and hence choose G 3 .In Equation 16, the h k 3,j are determined for the application.We seek to study the center manifold by reducing to an autonomous system and constructing an integral for the saddle such that, if set to zero, we obtain the center manifold.This means that we should choose h k n,j as follows.First, to eliminate time dependence, we choose the values h k n,j = 0 for j > 0 (j / ∈ J k ).Notice that by additionally setting h k n,0 = 0, we remove also the monomial associated with a k n,j .To construct the saddle integral, we remove all monomials with k 1 1 ̸ = k 2 1 .This means that the transformed Hamiltonian (up to some order N ) will have the form: where R is the remainder containing homogeneous polynomials of degree greater than N .Note that H depends on the product q 1 p 1 , rather than q 1 and p 1 separately.By defining the action variable I 1 and a symplectic conjugate hyperbolic angle θ 1 -see the Appendix in [32] for investigation of θ 1 -we obtain a canonical change of variables: If we neglect the remainder, the truncated Hamiltonian H N does not depend on θ 1 .Hence, I 1 is a first integral of the system, called the "saddle integral" [32].If we set I 1 = 0 (in particular, q 1 = p 1 = 0), then H N provides an N th -order approximation of the center manifold.We will exploit the saddle integral further in the next section to study the normal hyperbolic behavior of the center manifold.Note that we could additionally remove monomials to construct two integrals for the center manifold-this is the Birkhoff normal form.The benefit of the Birkhoff normal form is that it is integrable; however, the downside is the decreased radius of convergence.We note the construction of the center manifold is not unique.Depending on which monomials survive and which monomials are removed, one can obtain a different representation, e.g., the center manifold of L 2 in the QBCP is computed differently in [31] and [29].
Following this procedure, the variables of H N (0, q, p) will be complex variables in general.Hence, a canonical realification transformation is performed using the inverse of the complexification transformation given by Equation 10.Further, the nonlinear change of variables is obtained by the sequence of canonical transformations generated by G 3 , G 4 , . ... We perform this order-by-order and can transform each coordinate q j or p j via the transformation defined by where L n G k is the n-th order Lie derivative.See [33] for details.Finally, we investigate the dynamics in the local center manifold around DE L 4 by studying the vector field induced by H N (0, q, p) : R 4 → R which is independent of time.In these dynamics, the energy integral is recovered.This reduces the problem of studying dynamics around DE L 4 to the study of area-preserving maps parameterized by the energy [18].Once the energy is fixed, we can take a particular Poincaré section to obtain a 2-dimensional picture of the local center manifold.In this work, we use the vertical section Σ v = {q 2 = 0}, which corresponds to fixing z = 0 in synodical coordinates to first order [28].

Computation of Invariant Tori
In this work, we compute families of 2-dimensional quasi-periodic invariant tori by computing a parameterization v : T 2 → T ⊂ R 6 that maps angles θ = (θ 0 , θ 1 ) to a state on the quasi-periodic orbit, T , that is invariant under the flow induced by the Hamiltonian function H.The dynamics on the standard 2-torus, T 2 , are linear and the constant frequencies of oscillation, ω = (ω 0 , ω 1 ), are non-resonant.In particular, the frequency ω 0 is the frequency of the perturbation, i.e., ω 0 = ω S = 2.
Since the parameterization v is invariant under the flow, the invariant curve defined by the map u(θ 1 ) = v(0, θ 1 ) : where φ t is the flow induced by H for time t ∈ R, T = 2π/ω 0 = π is the period of the perturbation, and ρ = 2πω 1 /ω 0 is the rotation number of the torus.
To compute u explicitly, we discretize the invariant curve u(θ 1 ) over 2N + 1 angles whence the computation of invariant curves becomes a boundary value problem in which the 2N + 1 states must meet the boundary condition described by Equation 21.We solve the boundary value problem by defining the grid of angles, the grid of points along the invariant curve, and the free variable vector as respectively.Defining the constraint vector as where we apply the rotation operator R −ρ , sending arguments of functions depending on θ 1 to θ 1 − ρ, and the flow φ T is applied to each block element, as in [21], solutions to g(z) = 0 implicitly define a smooth manifold that is the 2-dimensional quasi-periodic invariant torus.It is important to note that the rotation operator can be constructed by composing discrete Fourier transform operations [21,27].Note that we solve this boundary value problem using a multiple shooting scheme with a Newton update step.
As any phase shift also satisfies the boundary value problem, we also add a phase constraint that minimizes the phase difference between the initial curve and the previously computed invariant curve, following [21,27].To compute families of invariant 2-tori, which lie in Cantorian oneparameter families in the HR4BP [34], we also employ a pseudo-arclength continuation scheme [23].Hence, we add the continuation constraint where Ũ and ρ denote the previously computed grid points along the invariant and rotation number, repsectively, n u and n ρ are approximations of the null space direction of the constraint Jacobian matrix corresponding to the invariant curve and rotation number, and ∆s is the step size.

Normal Behavior of Invariant Curves
We are interested in studying the stability of 2-dimensional quasi-periodic orbits in the HR4BP.By using the flow map method, we reduce the problem to that of computing the linearized normal behavior of an invariant curve of a diffeomorphism (the flow map φ T ) [35].Given an invariant curve u(θ 1 ) such that φ T (u(θ 1 )) = u(θ 1 + ρ) for all θ 1 ∈ T 1 .A small displacement h ∈ R 6 with respect to a point u(θ 1 ) on the curve is: The linear normal behavior is described by the system where A(θ 1 ) = D u φ T (u(θ 1 )).The rotation operator defined in the previous section is applied to the generalized eigenvalue problem as R ρ : ψ(θ) → ψ(θ +ρ).We consider the generalized eigenvalue problem in which we look for pairs (λ, ψ) To compute the stability of invariant tori, we discretize the operator R −ρ • A(θ 1 ), which is represented by a matrix.We then compute the eigenvalues and eigenvectors of this matrix by a standard numerical procedure.The resulting eigenvalues will form circles in the complex plane, covered by the same number of points as in the discretized invariant curve.If the invariant curve u(θ 1 ) is reducible, then its normal behavior will depend on six independent eigenvalues and eigenfunctions from this discretized set [35].As the map φ T is symplectic, the corresponding 6 eigenvalues will come in reciprocal and conjugate pairs, and there will be one pair of unity eigenvalues.The eigenfunctions corresponding to the unity pair are the tangent vector of the invariant curve along the flow and the tangent vector pointing in the direction of the one-parameter family of invariant 2-tori.
While the eigenvalues computed in this discretization will have different accuracy, due to the discretization, we implement a sorting method described in Jorba [35] and summarize the technique here.Representing each eigenfunction as a Fourier series j ψ j exp(ijθ 1 ), we want to determine which eigenfunctions have the least discretization error.Hence, we consider the truncated Fourier series The discretization error of the eigenfunctions is given by the magnitude of Fourier coefficients for j > N [21].As the eigenfunctions are analytic, their Fourier coefficients should decay exponentially from where the spectrum is centered [21].So, by checking the truncation error defined by we find the most "centered" eigenfunctions, and we choose these as our representative eigenfunctions and corresponding eigenvalues as the six representatives of the equivalence classes, provided the tails are sufficiently small.These representative eigenvalues and eigenfunctions then provide an accurate description of the normal dynamics.In particular, an eigenvalue pair on the unit circle in the complex plane describes a normally elliptic direction, and an eigenvalue pair on the real line away from ±1 describes a normally hyperbolic direction.Changes in stability type indicate a bifurcation in the family of quasi-periodic invariant tori, similar to Krein collisions of periodic orbits in Hamiltonian systems [36].

Hyperbolic Invariant Manifolds of Invariant Tori
The stable manifold theorem gives the existence of stable and unstable invariant manifolds emanating from quasi-periodic invariant 2-tori having a pair of real eigenvalues with modulus different from 1.The (global) stable and unstable manifolds are defined as the collections of points that asymptotically approach the invariant 2-torus as t → +∞ and t → −∞, respectively.In the section above, we described a process for computing eigenvalues and eigenfunctions of invariant curves.We follow [37,38] to compute the hyperbolic invariant manifolds of a normally hyperbolic invariant 2-torus.Suppose we have computed some invariant curve u(θ 1 ) which has a pair of hyperbolic eigenvalues λ s,u ∈ R with corresponding eigenfunctions ψ s,u (θ 1 ).For h ∈ R small and any θ ∈ [0, 2π], these hyperbolic eigenvalues and eigenfunctions satisfy: where A(θ 1 ) is defined in the previous section.Hence, the map is a parameterization of the linearization of the stable and unstable manifolds along the invariant curve u(θ 1 ).To compute the global stable and unstable manifolds of the invariant curve, one can fix a value of h = h 0 so that the error of the linearization is sufficiently small-we have used h 0 = 10 −5 so that the error is O(10 −10 ).However, as the HR4BP is π-periodic, under the π-periodic flow map a point of the invariant curve perturbed by h 0 onto the stable (unstable, respectively) manifold will be perturbed by h 0 /λ s (λ u h 0 , respectively).By choosing h ∈ [h 0 , h 0 /λ s ] ([h 0 , λ u h 0 ], respectively), we can effectively parameterize (in a linear sense) the stable and unstable manifolds of not only the invariant curve u(θ 1 ), but of the entire 2-torus v(θ 0 , θ 1 ) under the π-stroboscopic map.In other words, we have parameterized the fundamental cylinder of the hyperbolic invariant manifold.To globalize the manifold, we define a mesh of points along the cylinder (θ 1 , h) ∈ [0, 2π] × [h 0 , h 0 /λ s ] ([h 0 , λ u h 0 ], respectively) and propagate the points on the cylinder backward (forward, respectively) to compute the stable (unstable, respectively) manifold.Note that we perturb in ±h to obtain both sides of the invariant manifold.In this work, we choose a 251 × 251 grid for computations.

Results
In this section, we display and analyze the results from the methods described in the previous section applied to the Sun-Earth-Moon HR4BP.Firstly, we include a Poincaré map generated via center manifold reduction.Limitations of this method applied to L 4 are discussed.Then, five families of Lyapunov invariant 2-tori are presented along with their normal behavior.Finally, we show the existence of natural trajectories that escape the region around the triangular points, escaping the system or encountering regions near the primaries.

Poincaré Map
The center manifold reduction was computed about DE L 4 up to degree 12. Figure 4 shows the radius of validity for the center manifold reduction around DE L 4 up to degree 12. Observe that this radius decreases suddenly and asymptotically approaches zero.The small radius of convergence indicates that the Poincaré map will be severely limited in its utility to accurately describe dynamics near DE L 4 .Namely, Figure 5 shows a Poincaré map of the center manifold reduction of degree 4 at a fixed energy value h = 0.01.Note that the order and recovered energy integral are both small due to the divergence properties of the series expansion.Figure 5 shows the existence of 2-and 3-dimensional planar Lyapunov invariant tori around DE L 4 .Yet, the divergence properties of the series expansion limit the computation and numerical precision of the invariant tori.Hence, we utilize the flow map method to parameterize individual invariant 2-tori to study the invariant manifold structures around DE L 4 and L 4 2:1 resonant periodic orbit.

Families of Invariant 2-tori
In this section, we present the results of computing the five families of Lyapunov invariant 2tori around DE L 4 and the L 4 2:1 resonant periodic orbit and their normal behavior.We will identify bifurcations and compute bifurcations in the cases which are created by the perturbation.
In other words, if the bifurcation exists between periodic orbits in the Earth-Moon CR3BP, then the bifurcation will be identified but not computed or continued.

DE L 4 H
The first family of Lyapunov invariant 2-tori around DE L 4 is the planar family, which we denote as "DE L 4 H."The "H" stands for horizontal, as opposed to "V" for vertical.By taking as an initial guess the excited planar center mode of DE L 4 periodic orbit, we compute this family of invariant 2-tori using pseudo-arclength continuation and the algorithm described in Section 3.This family is qualitatively similar to (and a generalization of) the long-period family of periodic orbits around L 4 in the Earth-Moon CR3BP.
Figure 6 shows a hodograph of the family computation, along with sample torus representations.The family grows as x 0 increases.We follow the convention to represent families of invariant tori as with periodic orbits, i.e., to choose a particular state variable-x(τ ) at τ = 0 mod π-as a variable plotted against the frequency associated to the center mode.However, unlike in the case of periodic orbits, the choice of this variable is not unique-any point along the τ = 0 mod π invariant curve is equivalent.We choose the first x-coordinate along this invariant curve given by the numerical output of the differential correction.This results in a smooth curve for the family's hodograph.The family computation ends when the algorithm fails to converge.In this instance, the τ = 0 mod π invariant curve became less analytic, resulting in an increased number of Fourier modes used to describe the curve.While the family could be continued farther, we stopped the computation when  the maximum number of Fourier modes was reached.Alternatively, one could use the invariant curve along the other angular variable.
Figure 7 shows the computed normal behavior of the family of tori.As in the previous figure, the family grows in increasing x 0 = x(0 mod π).To show the stability, we describe two eigenvalues representing each normal mode.Due to the symplectic nature of the system, these eigenvalues come in reciprocal pairs (Krein quartets, in fact), so we pick one from each pair.The normal behavior for the DE L 4 H family remains of type center × saddle, as seen most clearly by the magnitude of eigenvalues λ 1 and λ 2 , the center and unstable eigenvalues, respectively.The behavior of this family is the most straightforward of the five presented, and we will return to DE L 4 H when considering transport near L 4 in the HR4BP in a later section.

DE L 4 V
The DE L 4 V family is a collection of invariant 2-tori emanating from the center mode of DE L 4 in the vertical (z) direction.This family is qualitatively similar to (and a generalization of) the vertical Lyapunov family of periodic orbits around L 4 in the Earth-Moon CR3BP.Using the same computational methods as above, we computed invariant 2-tori until a maximum z-value was reached because we are interested in the dynamics within this bound.
Figure 8 shows a hodograph of the family computation, along with sample torus representations.Note that the family grows as ω 2 increases.The invariant tori of this family retain a figure-8 shape as in the periodic vertical Lyapunov orbits about L 4 in the CR3BP; however, due to the periodic forcing of the Sun, an additional frequency is added.
Figure 9 shows the computed normal behavior of the family of tori.Unlike the DE L 4 H family, we show the normal behavior of DE L 4 V as a function of ω 2 because the family is one-to-one in this variable.Once again, note that the family grows as ω 2 grows.This figure shows that there are two bifurcations detected along the orbit family: a frequency-halving bifurcation, and a tangent bifurcation.The tangent bifurcation is seen in the CR3BP within the vertical Lyapunov family of periodic orbits around L 4 , and so we say that it is an artifact of the CR3BP and leave its computation for future work, if desired.The frequency-halving bifurcation is not an artifact of the CR3BP and is thus born from the periodic forcing of the Sun.The term "frequency-halving" is a generalization of "period-doubling" from periodic orbits to quasi-periodic orbits (which have no period but instead have multiple frequencies).We note that the frequency-halving bifurcation is detected because a Krein collision of eigenvalues occurs along the negative real axis.
We can summarize the normal behavior as follows.The family begins with the stability type of DE L 4 : saddle × center.There is a frequency-halving bifurcation that changes the stability of DE L 4 V to type center × center.Then, there is a tangent bifurcation which changes the stability to center × saddle.

L 4 2:1 H1
As the 2:1 resonant periodic orbit around L 4 is elliptic, there are two planar oscillatory modes.We adopt the naming convention of H1 and H2 for these horizontal families of invariant 2-tori.The L 4 2:1 H1 family is born from the planar oscillation corresponding to the "long-period" (i.e., smaller frequency) central mode.Thus, the H2 family is born similarly but from the "short-period" (i.e., larger frequency) central mode.
Figure 10 shows a hodograph of the family computation, along with sample torus representations.Note that the family grows as ω 2 decreases.Notice that the orbit appears to be getting pulled into DE L 4 .While it has been shown that the stable and unstable manifolds of DE L 4 do not intersect [10], i.e., DE L 4 is not a homoclinic point, one can apply the following co-dimension argument to see that we expect there to be homoclinic and heteroclinic 2-tori near DE L 4 .Consider an arbitrarily small planar 2-torus in the DE L 4 H family. Then the stable and unstable manifolds of this 2-torus intersect at a surface of dimension equal to the sum of the dimensions of the center-stable and center-unstable manifolds minus the dimension of the phase space.The center-stable and center-unstable manifolds each have dimension 3. Since we are considering the planar problem, the phase space is 5-dimensional (4 spatial, 1 temporal).Thus, 3 + 3 − 5 = 1, so we expect there to be a 1-parameter family of homoclinic connections in the planar problem.These homoclinic connections between invariant 2-tori in the plane act as partial barriers to transport.It is clear from the stroboscopic map in Figure 11 that the invariant curves are approaching a boundary close to the stable and unstable manifolds of DE L 4 .The authors speculate that the 2:1 H1 family ends in a homoclinic bifurcation with a homoclinic connection of an arbitrarily small member of the DE L 4 H family.The explicit computation of the connection is left for future work, as an identical argument shows 1-parameter families of heteroclinic connections exist, which may enable Arnold diffusion around DE L 4 in the Sun-Earth-Moon HR4BP.
Figure 12 shows the computed normal behavior of the family of tori.As in the previous figure, the family grows in decreasing ω 2 .Observe that this family is elliptic, i.e., has normal stability type center × center.Yet, as we continue along this family of elliptic tori, by analyzing the bottom-right plot showing arg(λ 2 ), we see that one pair of eigenvalues is approaching a Krein collision on the positive real axis, signaling an upcoming bifurcation.While we were unable to continue the family any further to exactly detect the bifurcation, this observation is consistent with the arguments made in the previous paragraph.

L 4 2:1 H2
As stated in the previous section, the L 4 2:1 H2 family of Lyapunov invariant 2-tori corresponds to the "short-period" (i.e., larger frequency) central mode of the elliptic 2:1 resonant periodic orbit.Figure 13 shows a hodograph of the family computation, along with sample torus representations.Note that the family grows as ω 2 decreases (or as x 0 increases).Unlike the 2:1 H1 family, the 2:1 H2 family is not pulled in by DE L 4 .Instead, the tori in this family appear to similarly wind around the 2:1 periodic orbit.We can see that the homoclinic and heteroclinic connections of the DE L 4 H family do not impose such transport restrictions on the H2 family as in the previous section with the H1 family.Figure 14 shows the invariant curves of H2 under the stroboscopic map.Note that as the colors move from blue to green to red, the family is increasing.Ultimately the computation of tori is stopped by the need to increase the number of Fourier modes above the maximum allowed in our computations-this threshold is reached due to the twisting of the invariant curves observed.
Figure 15 shows the computed normal behavior of the family of tori.As in the previous figure, the family grows in decreasing ω 2 .The family of invariant 2-tori inherits the normal behavior of the elliptic 2:1 resonant periodic orbit, i.e., the family begins with stability type center × center.There is one bifurcation identified in the family: a tangent bifurcation.This bifurcation connects the 2:1 H2 and V families, as can be seen in Figure 16 wherein one of the two broken branches bifurcates into the plane.The final of the five families of the Lyapunov 2-tori computed near L 4 is the 2:1 V family.Figure 16 shows a hodograph of this family with sample torus representations.First, observe that there is a symmetry from the left and right sides of the hodograph; because these tori emanate from a 2:1 periodic orbit, showing the τ = 0 mod π stroboscopic map yields two points representing the same torus on the hodograph.Moreover, the same principle applies to the previous two families, we simply omit one side because the two sides never meet, as in the 2:1 V family shown here.The family is computed from light to dark red.Due to the symmetry, the light blue and light green are equivalent representations of the same invariant torus (similarly from dark blue and dark green).The computation is stopped once the 2:1 periodic orbit is reached again.
Figure 16 shows a bifurcation in the 2:1 V family, and the blue and green torus representations indicate that these bifurcations connect the planar 2:1 H2 family, as seen in the previous section.It is important to note that this bifurcation is broken and that the light and dark blues thus represent different invariant tori (similarly with the light and dark greens, by symmetry).The precise cause of this bifurcation is unknown to the authors.The normal behavior was computed to be elliptic for all 2-tori computed, including the tori born from the bifurcation.The stability is not presented because, as the tori were all found to be elliptic, the equivalence class representatives were not chosen-no clear additional insight would be gained by their inclusion.
The orbit in the 2:1 V family which minimizes ω 2 (the torus represented in the middle subplot of Figure 16) is qualitatively similar to the DE L 4 V tori nearby-they are both in the shape of a figure 8.The frequency-halving bifurcation that occurs in the DE L 4 V family potentially relates the DE L 4 V and 2:1 V families of invariant tori.

Transport near L 4 in the HR4BP
Two of the central questions about the dynamics around EM L 4 in the Sun-Earth-Moon system are: can a point start close to DE L 4 and escape the system?How are trajectory fates dispersed?In this section, we show that the hyperbolic invariant manifolds of an invariant torus near EM L 4 can flow close to Earth and Moon, as well as escape the system.As this is qualitatively different behavior from the Earth-Moon CR3BP, we choose an quasi-periodic orbit near DE L 4 .In particular, we select a member of DE L 4 H, as shown in Figure 17.We take as the initial invariant curve to be at τ = 0 mod π, and we use h 0 = 10 −5 for the initial perturbation onto the stable or unstable manifolds.Figure 18 shows the results of cylinder set computations around an invariant torus of DE L 4 H close to DE L 4 .Blue points correspond to trajectories that escape the system; red points      intersect the radius of GEO in the configuration space; yellow points impact the lunar surface; black points do not satisfy any of the above criteria.The statistics for the cylinder sets are given in the table below:  The chaotic distribution of fates on each cylinder set can be described via closest approach to the Moon.We take as an illustrative example in the W u − map, choosing three points within a small neighborhood with each a distinct fate.Figure 19 shows the resulting trajectories (top) with the corresponding points on the cylinder (bottom).Note that the trajectories remain close together until the first lunar perilune, whence the trajectories enter the interior region near the Earth.Now, the yellow trajectory intersects the lunar surface.The red trajectory spends more time in the interior region then has a second close approach to the Moon and temporarily orbits the system before intersecting the GEO radius.Finally, the blue trajectory similarly winds around Earth in the interior region before making a close approach with the Moon-notably closer than the red trajectory-and picks up significant energy from the flyby causing it to escape the system.The change in energy from each lunar flyby can be seen by computing the Hamiltonian along each trajectory, as shown in Figures 20 and 21.The Hamiltonian is not constant along the flow as the system is periodic.Figure 20 shows the difference in energy change across each flyby, especially red and blue, i.e., GEO intersection and escapes the system.In the zoomed plot in Figure 21, we can more clearly see the effect of the first lunar perilune on the energy of each trajectory.Namely, the difference in energy along each trajectory is matched until the first lunar perilune, whence the energies are dephased and diverge from each other.Clearly, the yellow trajectory diverges least, as it intersects the lunar surface and hence has no lunar flyby.
Finally, as it is possible to move around the Earth-Moon system via stable and unstable manifolds of a quasi-periodic orbit close to EM L 4 , we can consider L 4 as a region of interest for the disposal of objects near the Moon.As we have seen, the stable manifold of the selected quasiperiodic orbit intersects the lunar surface at many points.We can compute the necessary ∆V for a single-impulsive maneuver sending a particle from the lunar surface directly onto the stable manifold of the chosen orbit, accounting for the velocity gained from the rotation of the Moon.This gives a conservative estimate on disposal maneuver cost for spacecraft along periodic and quasiperiodic orbits near L 1 and L 2 , as there is less energy difference between these orbits to apply a maneuver compared with the lunar surface.Figure 22 shows the comparison of all points along the stable manifold, W s ± , in time of flight (in years) and ∆V (in km/s), as well as the trajectory corresponding to the minimum time of flight.In the top plot, the red star is the minimum time of flight, and the blue star is the minimum ∆V .In the bottom plot, the first segment, shown in black, is the portion of the trajectory transiting into the greater L 4 region, while the second segment, shown in red, remains in the greater L 4 region.While the time of flight to the specific quasi-periodic orbit is nearly 10 years, the transit time to the greater L 4 region is approximately 78 days.Figure 25 shows a comparison between the continuation of L 4 in the BCP and in the QBCP.Due to their similarity in bifurcations and normal behavior, the literature considers the BCP as a sufficient approximation to model the dynamics around L 4 using the BCP or QBCP.
One key difference between the HR4BP and QBCP is the period of the perturbation.On one hand, the BCP and QBCP are periodic with period 2π/ω S , where ω S ≈ 0.9252 is taken to be the frequency of the Sun.On the other hand, the HR4BP is periodic with period π due to the symmetry provided by the Hill approximation of the Sun's gravitational effect.Consequently, the periodic orbit replacing L 4 in the BCP and QBCP will have period 2π/ω S , rather than period π as in the HR4BP.Similarly, the 2:1 resonant periodic orbit in the CR3BP is a different orbit when considering the periodic forcing of the QBCP or the HR4BP.In the HR4BP the 2:1 resonant orbit has period 2π in the CR3BP; in the QBCP the qualitatively similar resonant orbit would be 1:1 resonant orbit with period 2π/ω S .This CR3BP periodic orbit was not able to be continued into the Sun-Earth-Moon BCP.Despite this, the behavior of periodic solutions around L 4 is qualitatively identical between the BCP and HR4BP.

Invariant Tori
In [14], the vertical families of invariant 2-tori around PO1, PO2, and PO3 are computed in the Sun-Earth-Moon BCP.The planar families of invariant 2-tori around L 4 in the BCP have not been presented in the literature, though [7] computes some planar families of Lyapunov 2-tori in the BCP with the perturbation of solar radiation pressure.
Figure 26 shows the continuation of L 4 vertical families of invariant 2-tori in the BCP for several values of ε [14].The horizontal axis shows the value of ż of the invariant curve when z = 0, and the vertical axis shows the rotation number.In the Sun-Earth-Moon BCP, i.e., when ε = 1, there are 3 vertical families of invariant 2-tori emanating off of PO1, PO2, and PO3-these are denoted by F1, F2, and F3.Observe that there is a broken pitchfork bifurcation between the three families of invariant tori, similar to the bifurcation of periodic orbits.In fact, one can see from the continuation of these families in ε that the bifurcation of periodic orbits is the genesis of the relationship between the families of tori.The relationship between vertical tori, especially  the relationship between PO2/PO3 and PO1, is mirrored in the HR4BP; moreover, while the bifurcation in DE L 4 V has not been computed, preliminary calculations at smaller values of m suggest there may be a similar link relating the DE L 4 V and L 4 2:1 V families.The difference in the HR4BP is that the underlying bifurcation is not so clearly generated by the bifurcation of periodic orbits as in the BCP.This is under current investigation.

Powered Transfers to L 4
Finally, as we have discussed transport phenomena around L 4 and the utility for spacecraft disposal in cislunar space, we put this work into the context of the limited related literature.While there are no works considering transfers near the Moon to L 4 , there is some work considering transfers from Earth to L 4 in the BCP.Hence, a topic of continued study is to compute transfers between Earth and L 4 , which could be used for a communications satellite relaying information between Earth and the lunar gateway around EM L 2 .
Tan et al. consider single impulsive transfers from Earth using the window of easily approach (WOEA) without exploiting invariant manifolds-neither hyperbolic invariant manifolds nor quasiperiodic invariant tori [16].While the dynamics in the Sun-perturbed model allow for the exploitation of natural pathways along stable manifolds of invariant tori around DE L 4 , Tan et al. fail to take advantage of these natural structures which can alleviate some of the fuel cost in ∆V .Conversely, Liang et al. exploit quasi-periodic orbits near L 3 as parking orbits to compute powered transfers from Earth to L 4 [17].In their work, stable manifolds of several L 3 planar invariant tori are computed via backward propagation to find intersections with parking orbits around the Earth, ranging from low parking orbits to geostationary orbits.Unstable manifolds of the same L 3 planar invariant tori are then computed to identify intersections with stable regions around L 4 at a given epoch.A codimension argument shows that there are expected to be heteroclinic connections between 2-tori around L 3 and L 4 , but this is not considered in their work.

Conclusions
In this paper we have investigated the motion of a small particle moving near the Earth-Moon triangular points using as our model the Hill restricted 4-body problem (HR4BP).The periodic forcing of the Sun changes the L 4 equilibrium point into an isolated π-periodic periodic orbit, called DE L 4 .This perturbation qualitatively changes the linear behavior around the triangular points from elliptic to partially hyperbolic.Additionally, we show the continuation of the 2:1 resonant periodic orbit, explicitly computing the bifurcation killing one of the two branches that continue into the HR4BP.We considered the semi-analytical method of center manifold reduction, which we found to be extremely limited in utility.Hence, we computed 5 families of invariant 2-tori using the flow map method around DE L 4 and the 2:1 resonant periodic orbit, as well as their stability.A bifurcation relating to he vertical families of invariant 2-tori was identified.Finally, we investigated transport phenomena near EM L 4 by computing cylinder sets of a particular small planar quasi-periodic orbit close to DE L 4 .We compared results to existing literature of similar work in the bicircular restricted 4-body problem.

Fig. 3 :
Fig. 3: Stable and unstable manifolds of Branches B and A (DE L 4 ) under the π-stroboscopic map for several values of m.Note that the map shows τ = 0 mod π.

Fig. 4 :
Fig. 4: Radius of validity for center manifold reduction around DE L 4 .

Fig. 6 :
Fig. 6: Hodograph and sample invariant tori of DE L 4 H family. See text for details.

Fig. 7 :
Fig. 7: Normal behavior of DE L 4 H family of invariant tori.

Fig. 8 :
Fig. 8: Hodograph and sample invariant tori of DE L 4 V family.See text for details.

Fig. 9 :
Fig. 9: Normal behavior of DE L 4 V family of invariant tori.

Fig. 11 :
Fig. 11: Invariant curves of the π-stroboscopic map with W s,u of DE L 4 .

Fig. 14 :
Fig. 14: Invariant curves of the π-stroboscopic map.Note that W s,u of DE L 4 are not absolute barriers of transport.

Fig. 16 :
Fig. 16: L 4 2:1 V family with broken bifurcations.Note that there is a symmetry due to the 2:1 underlying period.

Fig. 18 :
Fig. 18: Sample stable manifold from a quasi-periodic invariant torus of the DE L 4 H family crossing GEO.

Fig. 19 :
Fig. 19: Trajectory comparison between adjacent points of W u − map with different fates.

Fig. 20 :
Fig. 20: Energy comparison between adjacent map points with different fates.

Fig. 21 :
Fig. 21: Zoom of energy comparison between adjacent map points with different fates.

Fig. 22 :
Fig. 22: ∆V [km/s] and time of flight [yrs] for direct transfers from the lunar surface onto the stable manifold of the selected small-amplitude quasi-periodic orbit near L 4 , including minimum time of flight trajectory.See text for details.

Fig. 26 :
Fig. 26: Continuation of L 4 vertical families of invariant 2-tori in the BCP for several values of ϵ [14].The horizontal axis shows the value of ż of the invariant curve when z = 0, and the vertical axis shows the rotation number.

Table 1 :
Stability of low-dimensional solutions in CR3BP and HR4BP