# Resonant quasi-periodic near-rectilinear Halo orbits in the Elliptic-Circular Earth-Moon-Sun Problem

Daniel Villegas-Pinto, Nicola Baresi, Slim Locoche, Daniel Hestroffer

► To cite this version:

Daniel Villegas-Pinto, Nicola Baresi, Slim Locoche, Daniel Hestroffer. Resonant quasi-periodic near-rectilinear Halo orbits in the Elliptic-Circular Earth-Moon-Sun Problem. Advances in Space Research, 2023, 71 (1), pp.336-354. 10.1016/j.asr.2022.08.011. hal-04011123

HAL Id: hal-04011123

https://hal.science/hal-04011123v1

Submitted on 16 Sep 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

ELSEVIER

Available online at www.sciencedirect.com

ScienceDirect

Advances in Space Research xx (2022) xxx-xxx

ADVANCES IN
SPACE
RESEARCH

(a COSPAR publication)

www.elsevier.com/locate/asr

# Resonant quasi-periodic near-rectilinear Halo orbits in the Elliptic-Circular Earth-Moon-Sun Problem

Daniel Villegas-Pinto<sup>a,*</sup>, Nicola Baresi<sup>b</sup>, Slim Locoche<sup>c</sup>, Daniel Hestroffer<sup>d</sup>

Received 12 April 2022;

## Abstract

Motivated by the near-future re-exploration of the cislunar space, this paper investigates dynamical substitutes of the Earth-Moon's resonant Near-Rectilinear Halo Orbits (NRHOs) under the Elliptic-Circular Restricted Four-Body Problem formulation of the Earth-Moon-Sun system. This model considers that the Earth and Moon move in elliptical orbits about each other and that a third body, the Sun, moves in a circular orbit about the Earth-Moon barycenter. By making use of this higher-fidelity dynamical model, we are able to incorporate the Sun's influence and the Moon's eccentricity, two of the most significant perturbations of the cislunar environment. As a result of these perturbations, resonant periodic NRHOs of the Earth-Moon Circular Restricted Three-Body Problem (CR3BP) are hereby replaced by two-dimensional quasi-periodic tori that better represent the dynamical evolution of satellites near the vicinity of the Moon. We present the steps and algorithms needed to compute these dynamical structures in the Elliptic-Circular model and subsequently assess their utility for spacecraft missions. We focus on the planned orbit for the NASA-led Lunar Gateway mission, a 9:2 synodic resonant $L_2$ southern NRHO, as well as on the 4:1 synodic and 4:1 sidereal resonances, due to the proximity to the nominal orbit and their advantageous dynamical properties. We verify that the dynamical equivalents of these orbits preserve key dynamical attributes such as eclipse avoidance and near-linear stability. Furthermore, we find that the higher dimensionality of quasi-periodic solutions offers interesting alternatives to mission designers in terms of phasing maneuvers and low-altitude scientific observations.

© 2022 COSPAR. Published by Elsevier Ltd All rights reserved.

**Keywords:** Quasi-periodic tori; near-rectilinear Halo orbit; Lunar Gateway; synodic and sidereal resonant orbits; dynamical systems theory; Restricted four-body problem

## 1. Introduction

The space sector has seen a returned interest in the manned and unmanned exploration of the Moon. In addition to the planned robotic lander missions by different space agencies, NASA's Artemis program will send humans back to the Moon's surface (NASA, 2020). Furthermore, a new NASA-led space station will be deployed about the Moon as part of NASA's Artemis program. This station, named the Lunar Gateway, will be placed in a lunar orbit and serve as an outpost for both lunar and interplanetary missions.

<footnote><sup>a</sup>PhD Candidate at CNES / IMCCE, Observatoire de Paris, Univ. PSL, CNRS, Sorbonne Univ., Univ. Lille. 77 av. Denfert-Rochereau, 75014 Paris, France</footnote>

<footnote><sup>b</sup>Surrey Space Centre, University of Surrey, Guildford, GU2 7XH, UK</footnote>

<footnote><sup>c</sup>Airbus Defence and Space, 31 Rue des Cosmonautes, 31400 Toulouse, France</footnote>

<footnote><sup>d</sup>IMCCE, Observatoire de Paris, Univ. PSL, CNRS, Sorbonne Univ., Univ. Lille, 77 av. Denfert-Rochereau, 75014 Paris, France</footnote>

<footnote>*Corresponding author.</footnote>

<footnote>E-mails: daniel.villegas-pinto@observatoiredeparis.psl.eu (D. Villegas-Pinto); n.baresi@surrey.ac.uk (N. Baresi); slim.locoche@airbus.com (S. Locoche); daniel.hestroffer@observatoiredeparis.psl.eu (D. Hestroffer).</footnote>
18

19

20

21

22

23

24

25

26

27

28

It will be used as a science and observation hub, for life support and habitation of visiting astronauts, and communications; it will
play a significant role in future programs and missions.

The Gateway's planned orbit type, a sub-family of the Earth-Moon's L₂ Southern Halo Orbit family known as the Near-Rectilinear Halo Orbit (NRHO) sub-family (NASA, 2019; Howell & Breakwell, 1984), has been extensively studied due to its relevance in the frame of space exploration. As reported in (Davis et al., 2017), the NRHO sub-family experiences lower perilune altitudes and presents attractive near-stable conditions with respect to the rest of the Halo family. The latter translates to lower station-keeping costs. Specifically, the targeted orbit for the Lunar Gateway is a 9:2 resonant synodic NRHO of the L₂ southern family NASA (2019), meaning that a spacecraft along this orbit completes 9 revolutions for every two revolutions of the Earth-Moon-Sun system, where the synodic period of the Earth-Moon-Sun is about 29.5 days. The use of synodic resonant orbits has been shown to be advantageous in the past, particularly regarding eclipse avoidance, as demonstrated in (Zimovan et al., 2017; Davis et al., 2017; Chikazawa et al., 2021).

It is common to study cislunar trajectories in the simplified Earth-Moon dynamical model of the Circular Restricted Three-Body Problem (CR3BP), where the NRHO and Halo families were originally formulated (Howell & Breakwell, 1984). However, the CR3BP can be a poor approximation of the real dynamics. For this reason, studies have been performed where NRHOs initially computed in the CR3BP are continued into high-fidelity ephemeris models and further analyzed in these models, as is the case in (Davis et al., 2017; NASA, 2019). Although multiple techniques exist, this is typically done by applying small velocity corrections (discontinuities) at specific points along the trajectory in order to keep it in the proximity of the baseline CR3BP orbit. Although the analyses performed via these techniques are very valuable, one could argue that part of their dynamical information might be lost due to the discontinuities applied, even if very small in magnitude. Moreover, we have no guarantees that the trajectories to which the algorithms converge represent the actual dynamical equivalent of the CR3BP's periodic NRHOs.

In order to retain a higher degree of dynamical information and a better match with respect to the CR3BP's trajectories, we can make use of higher-fidelity dynamical models that allow us to compute the real dynamical substitutes of the CR3BP's NRHOs. We can do this by using a model such as the periodic Bicircular Restricted Four-Body Problem (BCR4BP) (Gómez et al., 2001; Jorba-Cuscó et al., 2018), which incorporates the Sun gravity - a major influence in the dynamics of the cislunar space - with the Sun moving in a circular orbit around the Earth-Moon barycenter. In an effort to better capture the dynamical information of the NRHO sub-family and its Halo parent-family, previous work has focused on computing their dynamical substitutes in this model. In Boudad et al. (2020), synodic resonant periodic NRHOs are numerically continued into and analyzed in the BCR4BP model. Due to their resonant nature, these orbits remain periodic but complete a number of revolutions such that their period in the BCR4BP is equal to an integer number of system periods consistent with the resonance. The authors find that while general geometric properties are approximately maintained, such as apolune and perilune radii and eclipse avoidance properties, others, such as linear stability, are not necessarily conserved, which validates the need for analyses of the cislunar space in this type of higher-fidelity models. The work presented in Rosales et al. (2021) analyzes quasi-periodic tori families near the dynamical substitute of the L₂ point in the Earth-Moon-Sun BCR4BP model. Similar results are presented in Castellà & Jorba (2000), who first analyzed the vertical quasi-periodic families of the triangular equilibrium points in the same model. In McCarthy & Howell (2022), the authors compute transfers from Earth vicinity to quasi-periodic Halo orbits of the BCR4BP via invariant manifolds, whereas in Scheuerle & Howell (2021), low-energy impulsive transfers from the Earth vicinity towards the 9:2 synodic resonant NRHO and towards a Low Lunar Orbit (LLO) are computed in the BCR4BP.

While work in the literature remains very relevant and useful in the design of trajectories around the cislunar space, for the
most part it does not consider a significant perturbation present in the real dynamics of the system: the eccentricity of the Moon's
orbit. Aiming to further improve the accuracy of the dynamical systems used to model the cislunar environment, we propose to
incorporate this perturbation into our dynamical model, thus generalizing the BCR4BP into what we refer to as the Elliptic-Circular
Restricted Four-Body Problem (ECR4BP). The ECR4BP (or Elliptic-Circular) problem models the Earth and the Moon in elliptical
orbits about their barycenter and considers that the Sun moves in a circular orbit about the Earth-Moon's barycenter. The model
has been previously shown to describe the motion of a particle near the Earth-Moon's triangular equilibrium points accurately when
compared to JPL ephemeris, as opposed to the CR3BP (Castellà & Jorba, 2003). This is proven to be particularly significant in terms
of the practical stability of solutions, which the authors show the CR3BP and even the BCR4BP can over-estimate. Note that, just as
the BCR4BP, the ECR4BP is not a coherent model. That is, the motion of their massive bodies is not a solution of the Three-Body
Problem since the Sun is considered to not affect the motion of the Earth and the Moon. Nevertheless, the model remains a good
approximation of the dynamics of a spacecraft subject to the gravitational attraction of all bodies, making it an attractive option for
spacecraft mission design. While coherent restricted four-body models such as the Quasi-Bicircular model exist (Andreu, 1998),
we do not consider them in this work. Other restricted four-body models exist in literature, such as the Hill Four-Body Problem
Scheeres (1998); Scheeres & Bellerose (2005) or the BiElliptic model, which includes the ellipticity of all the primaries Assyrian &
Pourtakdoust (2010); Chakraborty & Narayan (2019). By making use of the ECR4BP to analyze the dynamical equivalents of some
of the CR3BP's NRHOs, among them, that planned for the Lunar Gateway, we aim to provide a closer match between the models
used by mission designers and the real dynamics of the system, thus offering important insights into the cislunar environment and
the possible trajectories of future spacecraft.

We propose to analyze three types of resonant NRHOs in this model, the 9:2 and 4:1 synodic resonant NRHOs, and the 4:1 sidereal resonant NRHO. In addition to the planned Gateway orbit, the two 4:1 resonances are chosen due their relative proximity to the 9:2 synodic resonance and their consistent stability behavior during continuation to other models (Boudad et al., 2020; Davis et al., 2017). Furthermore, although the synodic resonances have been studied in the Bicircular model, the sidereal resonance has only been analyzed in an ephemeris model Davis et al. (2017). Neither has been analyzed in a non-ephemeris dynamical model that incorporates both the gravitational attraction of the Sun and the eccentricity of the Moon's orbit.

Finally, we note that due the quasi-periodic nature of the Elliptic-Circular model, differently from the Bicircular case, periodic
orbits do not exist. Instead, the substitutes of the CR3BP's resonant periodic orbits will take the form of two-dimensional quasi-
periodic tori. In order to compute these quasi-periodic tori we employ a Two-Point Boundary Value Problem (TPBVP) algorithm
known as GMOS Gómez & Mondelo (2001); Olikara (2016); Baresi et al. (2018), which iteratively corrects an initial guess via the
Newton method until convergence criteria is verified. This allows us to obtain an initial invariant circle that can be integrated and
interpolated to obtain any state on the surface of the quasi-periodic torus. These quasi-periodic tori can then be analyzed in terms of
stability, geometry, and other relevant mission parameters. We highlight that, to the best of our knowledge, this is the first attempt
at calculating and analyzing quasi-periodic trajectories that take into account both the perturbation of the Sun and the eccentricity
of the Moon's orbit in the context of trajectory and mission design.

## 2. The Elliptic-Circular Restricted Four-Body Problem

The ECR4BP can be seen as a generalization of the Bicircular Restricted Four-Body Problem (Gómez et al., 2001), the latter of
which considers that the primary and secondary move in circular orbits whereas the former considers they move in elliptical orbits.
In both models, the third body (in our case, the Sun) moves in a circular orbit about the barycenter of the other two bodies, with all
bodies moving in the same orbital plane. The Bicircular model represents a periodic system, characterized by its synodic frequency or period, i.e. the rate at which the relative positions of the Earth-Moon-Sun repeat, which is given by $\Omega_3 = n - n_3$, where n is the mean motion of the primary and secondary and $n_3$ is the third body's mean motion. The Elliptic-Circular model, however, represents a quasi-periodic system with two (incommensurate) frequencies, one equal to $\Omega_3$, same as for the Bicircular, and another being the rate of the true anomaly of the two primaries, which is influenced by their eccentricity. Both models use a co-rotating reference

93

94

95

96

<figure><img src="image_3.png" /><figcaption>Fig. 1: Schematic representation of (a) the Bicircular Restricted Four-Body Problem and of (b) the Elliptic-Circular Restricted Four-Body Problem, in the co-rotating frame centered on the Earth-Moon barycenter. The dashed lines represent the orbits of the primary and secondary as viewed from a pseudo-inertial frame at the moment where the co-rotating and inertial frames coincide.</figcaption></figure>

97

98

99

100

101

102

103

frame centered on the barycenter of the primary and such that the x-axis constantly points from the primary to the secondary, the z-axis points in the direction of their mutual orbit's angular momentum vector, and the y-axis completes the orthogonal frame. A representation of each model can be seen in Figure 1 in the co-rotating reference frame centered on the barycenter of the primary and secondary. The dashed lines represent the orbits of the primary and secondary as viewed from a pseudo-inertial frame at the moment where the co-rotating and inertial frames coincide to illustrate the differences in the orbits of the primaries between the two models.

104

105

106

107

Since the distance between the two primaries is not constant in the ECR4BP, nor is their angular velocity, we employ the same normalization scheme used in the Elliptical Restricted Three-Body Problem (ER3BP) (Szebehely, 1967). The length unit is set to be equal to the distance between the two primaries, which arises from the equations of the two-body problem and which depends on the value of the true anomaly:

108

$$l = \frac{a(1 - e^2)}{1 + e \cos \nu}, \qquad (1)$$

109

110

where a is the semi-major axis of the primary and secondary's orbit, e is their eccentricity, and $\nu$ is the true anomaly. In order to compute the time unit, we first compute the time derivative of the true anomaly:

111

$$\dot{\nu} = \frac{dv}{dt} = \frac{(GM)^{1/2} (1 + e \cos \nu)^2}{a^{3/2} (1 - e^2)^{3/2}}, \qquad (2)$$

112

113

where $M = m_1 + m_2$, $m_1$ and $m_2$ are the masses of the primary and secondary, respectively, and G is the gravitational constant. We can now use the distance between the two primaries, l, the time derivative of the true anomaly, $\dot{\nu}$, and the sum of the masses
of the two primaries, M, to define the units of length, time, and mass as [LU] = l, [TU] = 1/\dot{v}, [MU] = m₁+m₂, respectively. Note that because we are, just like in the ER3BP, using the rate of the true anomaly in the time normalization, our equations of motion will be written (or differentiated) with respect to the true anomaly.

Although the Elliptic-Circular model was, to the best of the authors' knowledge, first mentioned in (Castellà & Jorba, 2003), the equations of motion of the model were not presented in said reference. In order to derive these equations, it is best to start from the equations of the ER3BP (Szebehely, 1967) and introduce the terms relative to the Sun, which we can directly import from the equations of the Bicircular model (Gómez et al., 2001). Accordingly, we can write the dimensional acceleration due to the third body's influence as

$$ \ddot{\mathbf{r}}_3^* = -\frac{Gm_3^*}{r_3^{*3}}\mathbf{r}_3^* - \frac{Gm_3^*}{a_3^{*3}}\begin{bmatrix} x_3^* \\ y_3^* \\ z_3^* \end{bmatrix}, \quad (3) $$

where the superscript * represents dimensional values, $\mu_3^*$ is the third body's gravitational parameter, $a_3^*$ is the semi-major axis of the third body's orbit, $x_3^*, y_3^*, z_3^*$ are the position coordinates of the third body in the co-rotating frame, $\mathbf{r}^* = [x^*, y^*, z^*]^T$ is the position-vector of the spacecraft, and $r_3^* = \sqrt{(x^* - x_3^*)^2 + (y^* - y_3^*)^2 + (z^* - z_3^*)^2}$ is the distance between the spacecraft and the third body. Note that the Sun is assumed to be on the same plane as the Earth-Moon system, which means $z_3^* = 0$. Then, using the normalization scheme,

$$ \ddot{\mathbf{r}}_3^* = \frac{GM}{l^2} \left( -\frac{\mu_3}{r_3^3} \mathbf{r}_3 - \frac{\mu_3}{a_3^3} \begin{bmatrix} x_3 \\ y_3 \\ 0 \end{bmatrix} \right) = \frac{GM}{l^2} \mathbf{r}_3'', \quad (4) $$

where □' denotes differentiation with respect to the true anomaly and $\mathbf{r}_3''$ is the normalized third-body acceleration. The term $\frac{GM}{l^2}$ simplifies with the other terms coming from the ER3BP normalization yielding the normalized equations of motion of the ECR4BP in the pulsating co-rotating frame:

$$ x'' - 2y' = \frac{\partial \psi}{\partial x}, \quad (5a) $$

$$ y'' + 2x' = \frac{\partial \psi}{\partial y}, \quad (5b) $$

$$ z'' = \frac{\partial \psi}{\partial z}, \quad (5c) $$

with

$$ \psi = \frac{1}{1 + e \cos v} \left[ \frac{1}{2} (x^2 + y^2 - ez^2 \cos v) + \frac{1 - \mu}{r_1} + \frac{\mu}{r_2} + \frac{\mu_3}{r_3} - \frac{\mu_3}{a_3^3} (x_3x + y_3y + z_3z) \right], \quad (6) $$

where $r_1 = \sqrt{(x + \mu)^2 + y^2 + z^2}$ and $r_2 = \sqrt{(x - (1 - \mu))^2 + y^2 + z^2}$. We can also write the equations of motion as

$$ x'' - 2y' = \frac{\partial \tilde{\psi}}{\partial x}, \quad (7a) $$

$$ y'' + 2x' = \frac{\partial \tilde{\psi}}{\partial y}, \quad (7b) $$

$$ z'' + z = \frac{\partial \tilde{\psi}}{\partial z}, \quad (7c) $$

with

$$ \tilde{\psi} = \frac{1}{1 + e \cos v} \left[ \frac{1}{2} (x^2 + y^2 + z^2) + \frac{1 - \mu}{r_1} + \frac{\mu}{r_2} + \frac{\mu_3}{r_3} - \frac{\mu_3}{a_3^3} (x_3x + y_3y + z_3z) \right], \quad (8) $$
which can be beneficial when taking partial derivatives, as will be shown in Section 4.1. Note that when setting the eccentricity to zero the equations simplify to the Bicircular case and when setting $\mu_3$ to zero they simplify to the ER3BP model. Due to the normalization used, the primary and secondary are always on the x-axis at $-\mu$ and $1 - \mu$, respectively.

To compute the position of the third body in the pulsating frame we use

$$ \begin{cases} x_3 = a_3 \cos \sigma = \frac{a_3^*}{l} \cos \sigma, \\ y_3 = -a_3 \sin \sigma = -\frac{a_3^*}{l} \sin \sigma, \\ z_3 = 0, \end{cases} \qquad (9) $$

where $\sigma$ is the angular position of the third body in the pulsating frame. Starting from

$$ \sigma = (v - n_3 t) + \sigma_0, \qquad (10) $$

we can derive the rate of change of the Sun's angular position with respect to the true anomaly in the pulsating co-rotating coordinate frame as

$$ \sigma' = \frac{d\sigma}{dv} = 1 - n_3/\dot{v} = 1 - \frac{n_3 (1 - e^2)^{3/2}}{n (1 + e \cos v)^2} \qquad (11) $$

with $n = \sqrt{GM/a^3}$. Owing to Eq. (11), it is found that $\sigma$ would grow linearly with the true anomaly $v$ as expected from Eq. (10). To enforce periodicity in $\sigma$, we find it best to introduce two auxiliary variables, namely $c = \cos \sigma$ and $s = \sin \sigma$. This is done in order for numerical algorithms to identify solutions separated only by a $2\pi$ variation on $\sigma$ as the same solution. Although modulo functions could be used with the same purpose, our implementation allows for derivatives and subsequently dependent equations to behave correctly numerically. Then, obeying to the equations of motion:

$$ c' = -s \sigma', \qquad (12a) $$

$$ s' = c \sigma'. \qquad (12b) $$

The latter may be rewritten in vector form as

$$ \chi' = \sigma' J \chi = \left( 1 - \frac{n_3 (1 - e^2)^{3/2}}{n (1 + e \cos v)^2} \right) J \chi, \qquad (13) $$

with $\chi = [c, s]^T$ and $J = \begin{bmatrix} 0 & -1 \\ 1 & 0 \end{bmatrix}$.

By appending Eq. (13) to the equations of motion Eq. (7) and $\chi$ to the six-dimensional state-vector $[r^T, v^T]^T$, we can formulate an eight-dimensional system of first-order ordinary differential equations such that

$$ X' = f(X, v, \xi) = f(X, v + 2\pi, \xi), \qquad (14) $$

where $X = [r^T, v^T, \chi^T]^T$ and $\xi = [\mu, e, \mu_3]^T$ is a vector of system parameters. In the following, periodic and quasi-periodic orbits are generated by first selecting sidereal and synodic Halo orbits in the CR3BP model of the Earth-Moon system ($\mu = 0.01215$, $e = 0.0$, $\mu_3 = 0.0$). These periodic orbits are later continued in either the ER3BP ($e = 0.0549$, $\mu_3 = 0.0$) or BCR4BP ($e = 0.00$, $\mu_3 = 328900.55$) by means of homotopy continuation (Boudad et al., 2020; McCarthy & Howell, 2022), and finally substituted by two-dimensional tori in the Earth-Moon-Sun Elliptic-Circular model via the GMOS algorithm (Olikara, 2016; Baresi et al., 2018; Villegas-Pinto et al., 2021) and pseudo-arclength continuation (Seydel, 2009).
### 3. Quasi-periodic tori

Computing a quasi-periodic torus is equivalent to finding the diffeomorphism $u(\theta) : T^d \to \mathbb{R}^n$ that maps points on the surface of a $d$-dimensional torus to states in the $n$-dimensional phase space of our dynamical system. The frequencies $\omega = \dot{\theta}$ represent the incommensurate constant torus frequencies and uniquely identify the quasi-periodic solution among different quasi-periodic candidates. A periodic orbit can be seen as a one-dimensional torus ($d = 1$), where $\theta = \theta_0$ represents the longitudinal angle along the curve and $\omega_0 = 2\pi/T$ is equivalent to the orbit's mean motion, where $T$ is the orbit's period. When increasing the dimension of the torus by one we are adding an angle coordinate, in this case transforming the our dynamical object to a two-dimensional torus such as the one of Figure 2. Figure 2 also illustrates an equivalent torus in $\mathbb{R}^n$ and the diffeomorphism $u$, which transforms between the two coordinate spaces. Each torus frequency gives us the rotation of their respective angle as a function of time. Because the frequencies are hereby assumed to be incommensurate, a trajectory on the surface of a quasi-periodic torus will cover its surface completely as time tends to infinity. A full $\Delta\theta_0 = 2\pi$ rotation along the $\omega_0$ direction takes a time $T = 2\pi/\omega_0$; we can then verify that after $T$, the angle $\theta_1$ will have rotated by $\rho_1 = \omega_1 T = \omega_1 2\pi/\omega_0$, where $\rho = \{\rho_1, ..., \rho_{d-1}\} = \{\omega_1 T, ..., \omega_{d-1} T\}$ is usually referred to as the rotation vector. Since $\rho_1$ is not a rational number, the value of $\theta_1$ after any number of revolutions $\Delta\theta_0 = 2\pi$ will never repeat. When considering two-dimensional tori it is common to refer to $\rho$ as the rotation number instead of rotation vector, as $\rho = \rho_1 = \rho$.

<figure><img src="image_4.png" /><figcaption>Fig. 2: Example of a two-dimensional torus in $T^2$ and the diffeomorphism $u$ that transforms the torus to and back from $\mathbb{R}^n$. The initial invariant circle is represented in red at $\theta_0 = 0$.</figcaption></figure>

The GMOS algorithm codified by (Gómez & Mondelo, 2001; Olikara, 2016; Baresi et al., 2018; Villegas-Pinto et al., 2021) is a Two-Point Boundary Value Problem (TPBVP) that makes use of the rotation vector and allows us to find the diffeomorphism $u$ of a torus by integrating a discretized initial invariant circle $X_0$ at some initial $\theta_0$ - see Figure 2 - over some time $T$ and reverting the rotation $\rho$ over its other angles. Indeed, for a two-dimensional torus, we can write

$$ \varphi_T(\mathbf{u}(\theta_0, \theta_1)) = \mathbf{u}(\theta_0, \theta_1 + \rho), \qquad (15) $$

where $\varphi_T$ represents the flow after a time $T$ and where we have the boundary conditions $\mathbf{u}(\theta_i + 2\pi, \cdot) = \mathbf{u}(\theta_i, \cdot)$. Note that $\varphi_T(\mathbf{u}(\theta))$ is the image of our stroboscopic map $u$ under a rotation along $\theta_0$ corresponding to $\Delta\theta_0 = 2\pi = \omega_0 T$. As such, we can define a rotation operator $[R_{-\rho}]$ that undoes the rotation along the other angular directions such that $[R_{-\rho}] \mathbf{u}(\cdot, \theta_1 + \rho) = \mathbf{u}(\cdot, \theta_1)$. The rotation operator can be written as a diagonal matrix by means of the DFT, as demonstrated in (Villegas-Pinto et al., 2021). We can then write the invariance condition for the GMOS algorithm for a discretized initial invariant circle $X_0$ integrated along the flow of the vector field over a time $T$:

$$ [R_{-\rho}]\varphi_T(X_0) - X_0 = 0. \qquad (16) $$
With the goal of finding and studying the dynamical substitutes of NRHOs in the ECR4BP, there are a few possible approaches that can be followed. However, it is important to first understand what these dynamical substitutes will be. In a quasi-periodic system such as the ECR4BP, a CR3BP periodic orbit that is non-resonant with either of the ECR4BP's perturbations would be replaced by a three-dimensional quasi-periodic torus that incorporates the frequencies of the system's perturbations and the period of the underlying periodic orbit. However, families of three- (and higher) dimensional quasi-periodic tori are harder to compute due to increased computational costs and more complex multi-parameter continuation algorithms (Henderson, 2002). If we choose underlying periodic orbits that are resonant with one of the system's perturbations, the tori become two-dimensional instead of three-dimensional, as one of the torus frequencies will incorporate both the period of the orbit and one of the system's periodic perturbations. As such, the solution envisioned in this study is to choose underlying periodic orbits of the NRHO family that are resonant with one of these perturbations, that is, either synodic resonant orbits (with the Earth-Moon-Sun period) or sidereal resonant orbits (with the Earth-Moon period), which will then produce two-dimensional resonant quasi-periodic tori in the Elliptic-Circular model.

## 4. Numerical continuation

The outset of our analysis are the CR3BP's $p:q$ resonant orbits, where $p$ is the number of orbital periods and $q$ is the number of either synodic or sidereal months. By design, synodic orbits remain periodic in the Bicircular model, whereas sidereal orbits remain periodic in the ER3BP. The difference with respect to their CR3BP counterparts is that the time-dependent periodic solutions will envision $p$-revolutions instead of just one in order to equal the $q$ revolutions of either the Bicircular or ER3BP system, respectively. A numerical continuation algorithm is implemented to compute these synodic and sidereal resonant orbits in their respective systems, while varying their corresponding dynamical parameter ($\mu_3$ for the synodic case and $e$, the Moon's eccentricity, for the sidereal one). This type of numerical continuation, where we decrease the number of parameters that describe our family by forcing and continue the family along one or more parameters that are external to object being computed, is typically known as homotopy continuation.

Both the homotopy continuation and GMOS algorithms can be implemented as a Two-Point Boundary Value Problem (TPBVP) that can be solved using Newton's method. In both cases, we can describe the problem as needing to satisfy boundary conditions at different points in time as

$$G(X_0, \xi) = 0, \qquad (17)$$

where $G: \mathbb{R}^N \times \mathbb{R}^k \to \mathbb{R}^N$, $X_0$ is our collection of initial states and $\xi$ are continuation and/or system parameters. Note that Eq. (17) can represent multiple shooting algorithms as well as single shooting algorithms. In the case of the trajectories computed in this paper, we use a multiple-shooting implementation as the one described in (Betts, 2010). Continuity between adjacent nodes is maintained by the condition

$$m^i(X) = \varphi_{\tau}(X^i, v^i) - X^{i+1} = X_{\tau}^i - X_0^{i+1}, \qquad (18)$$

where $i = \{0, 1, ..., N_0-2\}$, $N_0$ is the number of multiple shooting nodes, $\varphi_{\tau}(X^i, v^i)$ is the flow integrated from the initial conditions $(X^i, v^i)$ up to some time $\tau = T/N_0$ where $T$ is the total time. The flow of all our solutions after some time $\tau$ corresponds to the collection of states $X_{\tau}$. Note that Eq. (18) does not enforce any constraint between the last and the first node. In fact, this constraint
is enforced by another equation, g, such that G will be formed as

$$G = \begin{bmatrix} m^0 \\ \vdots \\ m^{N_0-2} \\ g \end{bmatrix} = 0. \qquad (19)$$

Typically, two different continuation algorithms are used, one for obtaining periodic orbits via homotopy continuation, where g usually takes the same form of Eq. (18) as $g = X_{\tau}^{N_0-1} - X_0^0$, and a different algorithm for obtaining quasi-periodic tori. However, we've verified that the numerical continuation algorithm hereby used to compute quasi-periodic can also be used, successfully, to perform homotopy continuation on periodic orbits, specifically for sidereal and synodic resonant orbits, from the CR3BP to the ER3BP and from the CR3BP to the BCR4BP, respectively. As mentioned, we make use of the Two-Point Boundary Value algorithm developed by Goméz, Mondelo, Olikara, and Scheeres (Gómez & Mondelo, 2001; Olikara, 2016; Baresi et al., 2018), which is typically known as GMOS, originally developed to compute and continue families of quasi-periodic tori. The GMOS algorithm uses the invariance equation presented in Equation (16) as g, such that $g = [R_{-\rho}]X_{\tau}^{N_0} - X_0^0$. When computing resonant periodic orbits, the rotation operator $[R_{-\rho}]$ does not affect $X_1$ due to the commensurate nature of the torus' frequencies associated with resonant periodic orbits, which allows us to continue these solutions as well. Additionally, due to the periodicity of these resonant orbits, the invariant circle collapses into a single point, which means we can set $N_1$, the number of points along the initial invariant circle, as $N_1 = 1$. Note that this continuation procedure is also effective with non-resonant periodic orbits, as $X_1$ will be at the “origin” of an invariant circle, and as such the rotation will have no effect.

Both continuation routines need a reliable initial guess; for the initial homotopy continuation, these are obtained by initializing the algorithm with the CR3BP's periodic orbit repeated over $p$ periods, whereas for obtaining the quasi-periodic tori we use the resonant multi-revolution periodic orbit previously obtained in the time-periodic systems, either the ER3BP or Bicircular model, to continue along the Sun mass or the eccentricity, respectively.

Starting from a known solution or good enough initial guess $(\tilde{X}_0, \tilde{\xi})$ we have that $G(\tilde{X}_0, \tilde{\xi}) = 0$. Then, by the Implicit Function Theorem, if

* $G_X(\tilde{X}_0, \tilde{\xi}) \in \mathbb{R}^N \times \mathbb{R}^N$ is non-singular

* and $G$ and $G_X$ are smooth near $(\tilde{X}_0, \tilde{\xi})$,

then there exists a unique continuous family $(X_0(\kappa), \xi(\kappa))$, along the continuation parameter $\kappa$, which may be one of the parameters in $\xi$, such that

* $(X_0(\kappa_0), \xi(\kappa_0)) = (\tilde{X}_0, \tilde{\xi})$,

* $G(X_0(\kappa), \xi(\kappa)) = 0$ for $\kappa$ near $\kappa_0$.

In order to find such families of solutions, we need to append additional equations to our TPBVP. These equations are known as phase and parametrizing conditions. The number of equations to be appended depends on the dimension of the manifold we are trying to compute and the system in which they are being computed. When computing $d$-dimensional manifolds in autonomous systems, we need to append $2d$ additional equations to account for the fact that these objects live in $d$-parameter families and are characterized by $d$ internal frequencies (Jorba & Villanueva, 1997; Muñoz-Almaraz et al., 2003; Villegas-Pinto et al., 2021). However, due to the non-autonomous nature of the systems considered and the fact that we are fixing one of the frequencies in order to choose resonant
objects, all the frequencies associated with each manifold will be fixed along the continuation procedures. This means that for both types of continuation routines described before (to obtain resonant periodic orbits and resonant quasi-periodic orbits), the solutions will lie in one-parameter families where the parameter is either the eccentricity of the Earth-Moon or the Sun mass parameter. It is worth noting that not only has the dynamical system been recasted as pseudo-periodic system via the incorporation of the Sun angle variables as part of the augmented state-vector, but also we are performing numerical continuation to obtain a family of solutions along the continuation parameter of interest. For that reason, although both of the tori's frequencies are fixed and known a priori, we include a phase condition associated with one of the torus' angular directions, as the ones described in (Villegas-Pinto et al., 2021; Olikara, 2016; Baresi et al., 2018)

$$p_0(X) = \frac{1}{N} \langle X, \frac{\partial \tilde{X}}{\partial \theta_0} \rangle = \frac{1}{N} \sum_{N_0-1}^{N_0-1} \sum_{N_1-1}^{N_1-1} X^T \frac{\partial \tilde{X}}{\partial \theta_0}, \qquad (20a)$$

$$p_1(X) = \frac{1}{N} \langle X - \tilde{X}, \frac{\partial \tilde{X}}{\partial \theta_1} \rangle = \frac{1}{N} \sum_{N_0-1}^{N_0-1} \sum_{N_1-1}^{N_1-1} (X - \tilde{X})^T \frac{\partial \tilde{X}}{\partial \theta_1}, \qquad (20b)$$

where $N = N_0 \times N_1$ for a two-dimensional torus. Whenever we perform continuation on the Sun mass, we use Eq. (20b), and whenever we continue the Earth-Moon eccentricity, we use Eq. (20a). Additionally, because the ECR4BP uses an augmented system and state-vector, we must add two phase constraints to fix the initial angle of the Sun, $\sigma_0$, at either 0 or $\pi$ as

$$p_c(X) = c_0 - \cos(\sigma_0), \qquad (21)$$

$$p_s(X) = s_0 - \sin(\sigma_0). \qquad (22)$$

For the parametrizing equation we include the pseudo-arclength constraint (Seydel, 2009), which is defined as

$$q(z) = \langle z - \tilde{z}, \frac{\partial \tilde{z}}{\partial h} \rangle - \delta h = \frac{1}{N} \sum_{N=0}^{N-1} (X_0 - \tilde{X}_0)^T \frac{\partial \tilde{X}_0}{\partial h} + (\xi - \tilde{\xi}) \frac{\partial \tilde{\xi}}{\partial h} - \delta s = 0, \qquad (23)$$

where $z = [X, \xi]^T$, $h$ is the pseudo-arclength step parameter, $\partial \tilde{z}/\partial h$ is the family tangent computed at the known solution, and $\delta h$ is the continuation step size. When the parameter we are continuing along is close to the desired value, we can switch from the pseudo-arclength constraint to a simple forcing equation in order to match said value. Finally, as explained in (Villegas-Pinto et al., 2021), we could use unfolding parameters to make the algorithm more robust. The full system with all the additional equations is then

$$F(X, \xi) = \begin{bmatrix} G(X, \xi) \\ p(X, \xi) \\ q(X, \xi) \end{bmatrix} = 0, \qquad (24)$$

where $p$ are the phase conditions and $q$ is the parametrizing condition. The system is solved iteratively via Newton's method as $DF(z)\delta z = -F(z)$, where $DF$ is the Jacobian matrix, until the norm of $F$ is below a specified tolerance, set here to $10^{-10}$. After converging on a family member, in order to produce the guess for the next member along the family we must compute the family tangent $\partial z/\partial h$ at the known solution from the null space of the $DF$ matrix, excluding the pseudo-arclength equation $q$. Then, the next family member is initialized as

$$z^{m+1} = z^m + \frac{\partial z^m}{\partial h} \delta h. \qquad (25)$$

In order to compute the Jacobian of $F$, we need to know the partials $\partial F/\partial X_0$, $\partial F/\partial \xi$, and $\partial F/\partial e$, where $\epsilon$ is the homotopy Sun mass parameter, which is multiplied by $\mu_3$ in the equations of motion and varied from zero to one throughout the continuation. These partial derivatives are calculated via variational equations as explained in the following subsection.
## 4.1. Partial derivatives with respect to continuation parameters

While the partial of $\partial F/\partial X_0$ is explained in (Olikara & Scheeres, 2012), the computation of $\partial F/\partial \varepsilon$ and $\partial F/\partial e$ partials for a system like the ECR4BP are not mentioned. In order to compute these quantities, we make use of variational equations, integrating the quantities of interest appended to the state-vector of our numerical integrator. We hereby focus specifically on $\partial G/\partial \xi$, given that the partials of $p(X, \xi)$ and $q(X, \xi)$ follow from the derivations in (Olikara, 2016; Baresi et al., 2018). Note that

$$ \frac{\partial}{\partial \xi} ([R_{-\rho}] X_{\tau}) = [R_{-\rho}] \frac{\partial X_{\tau}}{\partial \xi}, \qquad (26) $$

since $\rho$ does not depend on any parameter. As for $\frac{\partial X_{\tau}}{\partial \xi}$, we must obtain $\Theta = \partial\varphi(X_0)/\partial\xi$, where $\varphi(X_0)$ is the flow of the vector field $f = f(X, v, \xi)$ at the final state. We can find $\Theta$ by making use of the expression

$$ \Theta' = A\Theta + B, \qquad (27) $$

where $A = \partial f/\partial X$ and $B = \partial f/\partial \xi$, with initial conditions $\Theta_0 = 0$. Noting that the homotopy mass parameter $\varepsilon$ simply multiplies the normalized Sun mass in the equations of motion, the partial $\partial f/\partial \varepsilon$ can be derived directly from Eqs. (7a)-(8) as

$$ B = \frac{\partial f}{\partial \varepsilon} = \frac{1}{1 + e \cos \nu} r_3'', \qquad (28) $$

where $\tilde{\psi}$ is defined in Eq. (8). Conversely, the partial $\partial f/\partial e$ can be written as

$$ \frac{\partial f}{\partial e} = - \frac{\cos \nu}{(1 + e \cos \nu)^2} \nabla \tilde{\psi} + \frac{1}{1 + e \cos \nu} \frac{\partial r_3''}{\partial e}, \qquad (29) $$

where the dependency of $r_3''$ on $e$ comes from the normalization (note that $a_3 = a_3^*/l$ and that $r_3$ depends on $a_3$). The derivation of $\frac{\partial r_3''}{\partial e}$ can be computed by simply writing out the normalized length components relative to the Sun acceleration as function of the dimensional values divided by the length unit, as in Eq. (9). Using $g_{r_3} = \frac{-\mu_3}{r_3^3} r_3$ and $g_{a_3} = \frac{-\mu_3}{a_3^3} a_3$, with $a_3 = [c, -s, 0]^T = [x_3, y_3, z_3]^T$, we use the chain rule to write

$$ \frac{\partial r_3''}{\partial e} = \frac{dg_{r_3} dr_3}{dr_3 da_3} \frac{da_3}{de} + \frac{dg_{a_3} da_3}{da_3 de}. \qquad (30) $$

Both $dg_{r_3}/dr_3$ and $dg_{a_3}/da_3$ can be obtained from the expression

$$ \frac{dg_w}{dw} = - \frac{\mu_3}{w^3} \left[ I_3 - 3ww^T \frac{1}{w^2} \right], \qquad (31) $$

where $w$ must be replaced with either $r_3$ or $a_3$ to obtain the respective partial. Additionally, $dr_3/d a_3 = -I_3$, and for the last partial we have

$$ \frac{da_3}{de} = a_3 \left( \frac{2e}{1 - e^2} + \frac{\cos \nu}{1 + e \cos \nu} \right). \qquad (32) $$

The only dependency that cannot be computed explicitly is that of $c$ and $s$ with respect to $e$. These partial derivatives must be integrated along the trajectory as

$$ \frac{dc'}{de} = -s \frac{d\sigma'}{de}, \qquad (33) $$

$$ \frac{ds'}{de} = c \frac{d\sigma'}{de}, \qquad (34) $$

where $d\sigma'/de$ is obtained directly from Eq. (11) as

$$ \frac{d\sigma}{de} = \frac{n_3}{n} \left( \frac{3e(1 - e^2)^{\frac{1}{2}}}{(1 + e \cos \nu)^2} + \frac{2 \cos \nu (1 - e^2)^{\frac{3}{2}}}{(1 + e \cos \nu)^3} \right). \qquad (35) $$
326

## 4.2. Torus frequencies

327

As previously mentioned, the torus frequencies can be calculated and known a priori. This is due to the fact that we are computing resonant quasi-periodic tori in a quasi-periodic system that has been recasted as pseudo-periodic. The fundamental torus frequency, $\omega_0$, is computed from the period of the underlying (multi-revolution) periodic orbits as

328

329

330

$$ \omega_0 = 2\pi/T, \qquad (36) $$

331

where $T$ is the (normalized) period, which is resonant with either the synodic or the sidereal period. Although the expression for $\omega_0$ remains the same for the two types of resonances, the physical meaning of their torus angles and the expression for the second torus frequency, $\omega_1$, changes between the synodic and sidereal resonant quasi-periodic tori of the Elliptic-Circular model. These differences will be important when generating the guesses for the initial invariant circles of the different resonant tori.

332

333

334

335

For synodic resonant tori, the second torus angle, $\theta_1$, will reflect the true anomaly of the Earth-Moon, due to the additional perturbation imposed by the introduction of their eccentricity. For this reason, $\omega_1$ will be the normalized rate of the true anomaly in the system, which means $\omega_1 = 1$ for all synodic resonances, and, consequently, $\rho = T$, where $\rho$ is the rotation number. The fundamental torus angle, $\theta_0$, will reflect the motion of the Sun, where a full rotation over $\theta_0$ corresponds to the integer number of synodic periods that define the resonance, or in other words, $T = qT_{syn}$.

336

337

338

339

340

On the other hand, for sidereal resonant tori, since $T$ is resonant with the sidereal period, $\omega_0$ will be $1/q$, e.g. $\omega_0 = 1$ for a 4:1 sidereal resonance but for a 3:2 resonance we would have $\omega_0 = 0.5$. This means that the fundamental torus angle, $\theta_0$, will reflect the true anomaly of the Earth-Moon and that for a 4:1 and a 3:2 sidereal resonances a full rotation over $\theta_0$ would correspond to $\Delta\nu = 2\pi$ and $\Delta\nu = 4\pi$, respectively. The second torus angle, $\theta_1$, will then reflect the Sun angle for the sidereal resonances. However, because the rate of the Sun angle is not constant in the normalized system (recall that the time unit is not constant), we have to integrate the variation of the Sun angle over a full $2\pi$ rotation of $\theta_0$ to obtain:

341

342

343

344

345

346

$$ \omega_1 = \rho/T = \frac{1}{T} \int_{0}^{2\pi} \sigma' d\theta_0. \qquad (37) $$

347

It follows that, for the sidereal resonant case, $\rho$ will be equivalent to the variation of the Sun angle over one rotation of $\theta_0$.

348

## 5. Dynamical solutions

349

Using the GMOS algorithm with the previously described multiple shooting implementation, we can produce the dynamical solutions equivalent to the CR3BP's synodic and sidereal resonant periodic orbits of interest in the different dynamical systems. The parameters and constants used can be seen in Table 1. We point out that due to the continuation approach followed in this work, this choice of parameters could be different, and even applications to different systems could be undertaken using the same approach and setup hereby presented.

350

351

352

353

354

Additionally, it is worth noting that while all the resonances studied in this work converge correctly in the full Elliptic-Circular model using the aforementioned algorithm, many of the other synodic and sidereal resonant orbits fail to do so, even after being computed in the intermediate models, i.e. in the Bicircular and ER3BP models. A possible explanation to this phenomenon might be found in Boudad et al. (2020), where it is referred that convergence issues were also experienced when transitioning the 3:1 synodic resonant NRHO to a high-fidelity ephemeris model by the authors in (Davis et al., 2017). The reason put forward relates to the numerical continuation of this orbit and that of the 5:1 synodic resonant NRHO between the CR3BP and Bicircular, where a bifurcation along the stability curves can be seen at different values of the homotopy continuation parameter (this is the case

355

356

357

358

359

360
<table><caption>Table 1: Physical constants and parameters used for the different models. GM refers to the gravitational parameter, R refers to the body's radius, L to the distance or semi-major axis between two bodies. Sources: Wieczorek et al. (2006); IAU Division Working Group (2022); Mamajek et al. (2015).</caption><thead><tr><th>Parameter</th><th>Value</th></tr></thead><tbody><tr><td>GM Earth [m<sup>3</sup> s<sup>-2</sup>]</td><td>3.9860044189 × 10<sup>14</sup></td></tr><tr><td>GM Moon [m<sup>3</sup> s<sup>-2</sup>]</td><td>4.902801076 × 10<sup>12</sup></td></tr><tr><td>GM Sun [m<sup>3</sup> s<sup>-2</sup>]</td><td>1.32712440018 × 10<sup>20</sup></td></tr><tr><td>R Earth [km]</td><td>6378.137</td></tr><tr><td>R Moon [km]</td><td>1737.103</td></tr><tr><td>R Sun [km]</td><td>6.9551 × 10<sup>5</sup></td></tr><tr><td>L Earth-Moon [km]</td><td>384 399</td></tr><tr><td>L Sun-Earth [km]</td><td>149.5978707 × 10<sup>6</sup></td></tr></tbody></table>

for many resonances but extensive testing was not performed). This causes a mismatch in the eigenstructure of the intermediate multi-revolution periodic orbits, which seems to later prevent the numerical continuation to the full Elliptic-Circular model as quasi-periodic tori or to ephemeris models. Unlike other resonances, the ones selected in this paper do not present such a bifurcation or mismatch in the continuation to the intermediate models.

## 5.1. Initial guess generation

As the initial guesses of the quasi-periodic tori continuation procedures, the multi-revolution resonant periodic orbits must first be obtained in the intermediate models, i.e. in the Bicircular and ER3BP models. To do so, we start from the CR3BP resonant *p*:q-type periodic orbits and initialize the continuation procedure by repeating the orbits over *p* periods of the original CR3BP trajectories, such that their period will be *T* = *p* *T*<sub>res</sub>, where *T*<sub>res</sub> is the period of the resonant periodic orbits in the CR3BP model. As mentioned, we also make use of the GMOS algorithm to continue these resonant periodic orbits into the intermediate models. To do so, we set *N*<sub>1</sub> = 1, and, as we are using a multiple shooting implementation, set *N*<sub>0</sub> = *p*, such that there is only one multiple shooting node per revolution and that each one is placed at the apolune of its orbital revolution. We set a phase constraint with respect to the *θ*<sub>0</sub> angular coordinate for the continuation along Sun mass parameter (CR3BP to BC4RBP) or with respect to *θ*<sub>1</sub> for the continuation along the eccentricity (CR3BP to ER3BP) as detailed in (Villegas-Pinto et al., 2021), the two additional phase constraints on the Sun angle, and the pseudo-arclength equation, just as explained in Section 4. The initial resonant periodic NRHOs in the CR3BP and the result of this continuation procedure can be seen in Figures 3 and 4, respectively, where the 4:1 and 9:2 synodic and 4:1 sidereal resonant NRHOs are shown first in the CR3BP and then in their respective intermediate models, the Bicircular model for the synodic resonances and the ER3BP for the sidereal resonance. The trajectories are shown in dimensional coordinates in the Earth-Moon co-rotating reference frame, albeit centered on the Moon. The projections of each trajectory on the different planes can be seen in grey.

It is important to point out that all the solutions hereby presented, including the quasi-periodic solutions to be shown, were computed and numerically continued by choosing initial conditions or multiple shooting nodes at apolune points. As shown in (Davis et al., 2017), NRHOs are highly sensitive to errors and perturbations at perilune, which seems to be caused by strong non-linearities present around this region (Jenson & Scheeres, 2021). In the same way that these authors advise against applying correction maneuvers around perilune due to the high-sensitivity of said region and instead recommend applying them at apolune, the iterative correction algorithms hereby employed use as initial conditions states that are situated at apolune, thus applying the corrections to these state-vectors at the least sensitive positions along the trajectories. We find that this significantly improves the
<figure><img src="image_5.png" /><figcaption>Fig. 3: The resonant periodic southern NRHOs under study, computed in the CR3BP model and shown in the Earth-Moon co-rotating frame centered on the Moon in dimensional coordinates: (a) the 4:1 synodic resonant NRHO, (b) the 9:2 synodic resonant NRHO, and (c) the 4:1 sidereal resonant NRHO.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 4: The multi-revolution resonant periodic southern NRHOs under study, computed in their intermediate models and shown in the Earth-Moon co-rotating frame centered on the Moon in dimensional coordinates: (a) the 4:1 synodic resonant NRHO in the BCR4BP, (b) the 9:2 synodic resonant NRHO in the BCR4BP, and (c) the 4:1 sidereal resonant NRHO in the ER3BP.</figcaption></figure>

convergence of the algorithms, which might otherwise fail.

Additionally, it is worth noting that the synodic resonant solutions hereby presented are computed for initial Sun angles equal to zero, i.e. when all bodies are aligned and the Sun is in the positive x-axis. Likewise, the sidereal resonant solutions are computed by setting the initial true anomaly to zero. Nevertheless, note that both intermediate periodic systems, the Bicircular and the ER3BP, present the following symmetries Oshima (2022); Campagnola et al. (2008):

$$ (x, y, z, \dot{x}, \dot{y}, \dot{z}, \tau) \rightarrow (x, -y, -z, -\dot{x}, \dot{y}, \dot{z}, -\tau), \qquad (38a) $$

$$ (x, y, z, \dot{x}, \dot{y}, \dot{z}, \tau) \rightarrow (x, -y, z, -\dot{x}, \dot{y}, -\dot{z}, -\tau), \qquad (38b) $$

where $\tau$ represents the time-periodic variable of the system, which for the Bicircular is the Sun angle $\sigma$ and for the ER3BP is the true anomaly $\nu$. This indicates that a sufficient condition for symmetric periodic orbits in these models is that they cross the x-axis or xz-plane perpendicularly at either $\tau = 0$ or $\tau = \pi$. When considering $p:q$ resonant orbits with odd $p$, the orbit will have completed a non integer number of revolutions, namely $p/2$, after a time $q\pi$. Indeed, the different initial angle of zero or $\pi$ produce two different periodic orbits, which is not the case for even $p$. Although in this work we focus only on solutions with an initial angle of zero, one could extend these to cover initial angles of $\pi$ as well. Moreover, depending on the $p:q$ resonant orbit being
computed, one could consider other initial angles separated by a phase of $q2\pi/p$. However, this would lead to the same periodic orbit, simply initialized at a different point along its trajectory. An example would be the 4:1 multi-revolution synodic resonant NRHO in the Bicircular, which as shown in (Boudad et al., 2020), could be computed with initial Sun angles separated by a phase of approximately $360^\circ/4 = 90^\circ$, which all correspond to the same orbit. Previous authors have presented unique solutions in the Bicircular for different initial Sun angles, but this typically entails small discontinuities in the trajectories (Scheuerle & Howell, 2021). Since we admit a solution as converged when the error vector is below a certain threshold, using the Newton method instead of, for instance, least-squares, our algorithm does not accept these solutions, causing the error vector to plateau and not converge. Nevertheless, as mentioned, the solution space of the trajectories presented in this study could be extended by changing the initial angles of the intermediate periodic models to $\pi$, or even to include the discontinuities associated with different angles.

## 5.2. Quasi-Periodic Tori Continuation

After obtaining the multi-revolutions periodic solutions in the intermediate models, we are in a position to proceed with the numerical continuation towards the Elliptic-Circular model taking these solutions as initial guesses of the quasi-periodic tori. We take the converged solutions of the intermediate models, which are already defined in $N_0 = p$ multiple shooting nodes situated at the apolune of each revolution, and we form the invariant circles at each $\theta_0$ node by repeating each of these states $N_1$ times. Of course, this means that for the zeroth family member of this numerical continuation, i.e. where $e = 0$ for the synodic resonances or $\epsilon = 0$ for the sidereal resonance, all cartesian coordinates of the initial invariant circle will coincide with initial state of the underlying multi-revolution resonant periodic orbit. However, as we move along the family they will naturally spread and display the closed curve more clearly. We set the value of $N_1$ to 50 for higher accuracy but note that it can be set between 30 and 50 (or larger), depending on the resonant trajectories. We find that trajectories with lower perilunes, such as the 9:2 synodic resonant NRHO, typically require a higher number of $N_1$ nodes to achieve good quadratic convergence and to accurately interpolate points on the surface of the torus, although this is not strictly necessary to find the solutions. This also falls in line with the previous observations about the strong non-linearities present at perilune.

It should be noted that for the initial guess of the first family member, although each multiple shooting node will be comprised of invariant circles that share the same initial position and velocity vectors (since we simply repeated each state $N_1$ times for each node), the integration of their trajectories will not be equal. In fact, since we are performing continuation along either the Earth-Moon eccentricity for synodic resonant tori or Sun mass parameter for sidereal resonant tori, the $\theta_1$ angular direction of the torus will reflect different quantities, just as mentioned in Section 4.2. Specifically, it will reflect the Earth-Moon true anomaly or the Sun angle, respectively.

For a synodic resonant torus, each state along the initial invariant circle will have a different value of true anomaly, corresponding to the different values of $\theta_1$ as

$$v^{0,j} = 2\pi j/N_1, \qquad (39)$$

where $j = \{0, 1, ..., N_1 - 1\}$ and $v^{0,j}$ represents the $j$th of value true anomaly, or, equivalently, of $\theta_1$, of the initial node $i = 0$ along $\theta_0$. In the same way, the values of true anomaly for the subsequent nodes will be

$$v^{i,j} = Ti/N_0 + 2\pi j/N_1, \qquad (40)$$

with $i = \{0, 1, ..., N_0 - 1\}$.
436

437

For the initial guess of the first family member of the sidereal resonant tori, we vary $\theta_1$ along with the Sun angle $\sigma$, meaning that, for the initial guess, each state along the initial invariant circle will have a different Sun angle

438

$$ \sigma^{0,j} = 2\pi j / N_1. \qquad (41) $$

439

440

441

442

However, it should be noted that for sidereal resonant tori we cannot compute the initial guesses for the values of $\sigma^{i,j}$ at the other $N_0 - 1$ nodes in the same way as we do for $\nu^{i,j}$ in the case of synodic resonant tori. The reason for this is that $\sigma$ does not vary linearly with the true anomaly. Instead, we must first integrate $\sigma'$ over $\Delta\theta_0 = 2\pi$ and evaluate $\sigma$ at the $N_0$ nodes, similarly to how we obtain $\omega_1$ and the rotation number in Eq. (37). The values of $\sigma$ at the first point along each invariant circle will then be

443

$$ \sigma^{i,0} = \int_{0}^{2\pi i / N_0} \sigma' \bigg|_{\sigma=\sigma^{0,0}} d\theta_0, \qquad (42) $$

444

where $\sigma'$ is given by Eq. (11). Then,

445

$$ \sigma^{i,j} = \sigma^{i,0} + \sigma^{0,j}. \qquad (43) $$

446

447

448

449

450

Because $c$ and $s$ have been appended to the state-vector of our equations of motion, the values of $\sigma^{i,j}$ are left to vary freely and are corrected at each iteration of the continuation algorithm, constrained only by the multiple shooting equations that enforce continuity at the nodes and the constraint on the initial Sun angle of the initial invariant circle (cf. Eqs. (21) and (22)). As such, the expression on Eq. (43) is used only as an initial guess. Additionally, although $\sigma$ relates to $\theta_1$ in the case of sidereal tori, we see that they do not necessarily equal each other, as $\theta_1$ varies linearly with $\theta_0$, which is not the case for $\sigma$.

451

452

453

454

455

456

457

458

459

460

461

462

463

464

465

466

467

468

469

By considering the solution where the homotopy continuation parameter is zero as the first known solution, i.e., by repeating the trajectories coming from the intermediate models over $N_1$ points, we can take a step along the family tangent direction as described in Eq. (25). This is achieved by first computing the family tangent from the null space of the reduced DF matrix, i.e., the Jacobian matrix without the pseudo-arclength continuation equation. Then, by means of the GMOS algorithm and the proposed continuation procedure, we continue the solutions to produce resonant two-dimensional quasi-periodic tori up to an Earth-Moon eccentricity value of $e = 0.0549$ and a Sun mass homotopy parameter of $\varepsilon = 1$, corresponding to the actual Sun mass. The results of this continuation procedures can be seen in Figures 5 to 7, where the initial guesses used, i.e. the multi-revolution periodic orbits from the intermediate models, can be seen on the left and the quasi-periodic tori on the right. The solutions are plotted in dimensional coordinates in the Earth-Moon co-rotating frame centered on the Moon, and their projections on the different planes are visible in grey. Additionally, we can also look at the initial invariant circles of the quasi-periodic tori throughout the continuation procedure. This can be seen in Figures 8 and 9, which display said invariant circles in normalized and dimensional coordinates, respectively. As mentioned earlier, the nodes of the initial invariant circles progressively move away from the point that represents the resonant periodic orbits of the time-periodic intermediate systems as we move along the numerical continuation. Moreover, it is interesting to note that particularly for the case of the 4:1 synodic resonance, the dimensionalization of the position coordinates plotted with respect to the Moon, i.e. where the Moon is at the origin of the frame, deforms the invariant circles of Figure 9b when compared to those in normalized coordinates from Figure 8b. This is expected due to the previously mentioned values of true anomaly for the different nodes of the invariant circles of the synodic resonant tori (see Eqs. (39), (40)), which affect the value of the length constant, as seen in Eq. (1). The deformation is mostly seen along the x-axis, due to varying Earth-Moon distance for different values of true anomaly.
<figure><img src="image_7.png" /><figcaption>Fig. 5: (a) The periodic resonant synodic 9:2 NRHO in the Bicircular model and (b) the corresponding quasi-periodic tori of the 9:2 synodic resonant NRHO in the full ECR4BP model for a value of eccentricity *e* = 0.0549.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 6: (a) The multi-revolution periodic resonant synodic 4:1 NRHO in the BCR4BP model and (b) the corresponding quasi-periodic tori of the 4:1 synodic resonant NRHO in the full ECR4BP model for a value of eccentricity *e* = 0.0549.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Fig. 7: (a) The periodic resonant sidereal 4:1 NRHO in the ER3BP model for a value of eccentricity *e* = 0.0549 and (b) the corresponding quasi-periodic tori of the 4:1 sidereal resonant NRHO in the full ECR4BP model for the same value of eccentricity.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 8: The initial invariant circles of the resonant quasi-periodic tori for different family members along the numerical continuation in the ECR4BP, shown in normalized coordinates in co-rotating reference frame centered on the Moon. Initial invariant circles of the (a) 9:2 synodic resonant quasi-periodic NRHO (b) 4:1 synodic resonant quasi-periodic NRHO (c) 4:1 sidereal resonant quasi-periodic NRHO.</figcaption></figure>

Looking at the quasi-periodic tori, we can see that adding the second perturbation to the system translates into a “thickening” of the solution around the periodic resonant NRHO. Interestingly, we see that this thickening is most prevalent in the 9:2 synodic and 4:1 sidereal resonant quasi-periodic NRHOs. Additionally, when visually comparing the 9:2 synodic resonant NRHO to the long-term high-fidelity ephemeris model converged trajectories from literature (NASA, 2019; Davis et al., 2017), the quasi-periodic torus shown in Figure 5b seems thicker. Of course, even if both solutions have the same starting point – the NRHO of the CR3BP –, they are computed via different methods. The full-ephemeris solutions are integrated backwards and forwards in time for several years (10 or more) and present velocity discontinuities where correction maneuvers are executed. Even so, the natural quasi-periodic solution hereby presented is integrated for a much shorter period (less than a year), and it already displays a wider and seemingly thicker visual appearance. This suggests that the quasi-periodic trajectories hereby computed and the converged trajectories obtained
<figure><img src="image_11.png" /><figcaption>Fig. 9: The initial invariant circles of quasi-periodic tori for different members along the numerical continuation in the ECR4BP, shown in dimensional coordinates in co-rotating reference frame centered on the Moon. Initial invariant circles of the: (a) 9:2 synodic resonant quasi-periodic NRHO (b) 4:1 synodic resonant quasi-periodic NRHO (c) 4:1 sidereal resonant quasi-periodic NRHO.</figcaption></figure>

in full-ephemeris models might represent different dynamical solutions, although analyses of these quasi-periodic solutions in a full-ephemeris model would be necessary for a more accurate comparison.

## 6. Stability analysis

In order to better analyze the two-dimensional quasi-periodic tori obtained in the full Elliptic-Circular model, we study the linear stability associated with these dynamical solutions. Although different criteria exist in order to study the stability of dynamical solutions, we choose to analyze the solutions hereby presented using finite Lyapunov exponents. The reason behind this choice of metric is that it normalizes the stability values by the period of the solution. This is convenient due to the long periods introduced with the multi-revolution periodic NRHOs computed in the intermediate models. If another approach was used, such as the norm of the eigenvalues of the state transition matrix, the longer periods of these orbits would make it seem like their stability had significantly decreased with respect to their single-revolution counterparts of the CR3BP. By using Lyapunov exponent we obtain a more accurate interpretation of the changes in stability, as demonstrated by Boudad et al. (2020) when comparing the stability of the CR3BP resonant NRHOs with the synodic resonant multi-revolutions NRHOs of the Bicircular model. The Lyapunov exponents can be obtained as

$$ \phi_i = \mathcal{R} \left( \frac{\ln \lambda_i}{T} \right), \qquad (44) $$

where $\lambda_i$ are the eigenvalues of the state transition matrix after one period, i.e. the monodromy matrix, and $T$ is the period of the solution. A solution is considered to be linearly stable if all its Lyapunov exponents are equal to zero. It should be noted that this criteria is formulated for periodic orbits, and that it arises from a limit approximation (Scheeres, 2012). Although we do not derive said approximation to employ the Lyapunov exponents in the stability analysis of quasi-periodic solutions, we tacitly use this metric and find that, empirically, it produces consistent results that are congruous with those shown for their underlying periodic orbits. As such, we obtain the eigenvalues associated with the Floquet matrix, $B$, of the quasi-periodic tori and use the tori's fundamental period, i.e. one full revolution along $\theta_0$, which corresponds to the period of the underlying multi-revolution periodic NRHO, to compute their Lyapunov exponents. As explained in (Olikara, 2016; Baresi, 2017), the eigenvalues of a quasi-periodic torus' Floquet matrix provide us with linear stability information about the torus, similarly to the monodromy matrices of periodic
orbits. The Floquet matrix is computed by assembling a block diagonal of the $N_1$ monodromy matrices, $\Phi_i$, of each trajectory along the invariant circle, and undoing the torus rotation by means of the rotation operator, $[R_{-\rho}]$, as

504

$$ \mathbf{B} = [R_{-\rho}] \begin{bmatrix} \Phi_0 & & & \\ & \Phi_1 & & \\ & & \ddots & \\ & & & \Phi_{N_1-1} \end{bmatrix}, \qquad (45) $$

which will have a dimension of $nN_1 \times nN_1$, where n is the number of state variables. Since the stability of the Sun's position bears no physical meaning, we can focus only on the first six rows and columns of each monodromy matrix. We can then compute the eigenvalues of the Floquet matrix, which will be distributed in concentric circles on the complex plane. Since eigenvalues occur in reciprocal pairs, so will the radii of these circles, i.e. $R_1 = 1/R_2$, where $R_1$ and $R_2$ are the radii of reciprocal circles in the complex plane (Jorba & Villanueva, 1997). Due to the six-dimensional state we are considering for the stability analysis of these quasi-periodic tori, we will have six concentric circles, where each is a reciprocal pair of another. We can then either sample each circle to obtain an eigenvalue or simply consider their radius, which will be equal to its intersection with the real axis, in order to obtain the six $\lambda_i$ necessary for the computation of the Lyapunov exponents.

As mentioned in the beginning of Section 5, Boudad et al. (2020) perform the same analysis when translating periodic synodic resonant NRHOs from the CR3BP to the Bicircular model, and find that out of the four synodic resonances studied by the authors, the 3:1 and 5:1 synodic resonant NRHOs present bifurcations along the stability curve of their Lyapunov exponents when plotted against the continuation parameter. These bifurcations, which change the stability of these two NRHOs from linearly stable to slightly unstable, seem to translate in convergence issues when continuing the solutions to high-fidelity ephemeris models. We verify that the same resonant orbits also experience convergence issues when continuing them into the full Elliptic-Circular model from the Bicircular. Accordingly, the synodic resonances hereby considered (9:2 and 4:1), which as shown by Boudad et al. do not present any bifurcations along their stability curve when continuing them from the CR3BP to the Bicircular model, do not produce convergence issues in their continuation to the full Elliptic-Circular model.

Furthermore, the analysis of the Lyapunov exponents of quasi-periodic solutions as a function of their continuation parameter allows us to verify that the solutions correspond to the actual dynamical equivalents of the CR3BP's solutions by ensuring that there are no discontinuities along the family curve.

The evolution of the Lyapunov exponents of each of the resonant quasi-periodic NRHOs along their continuation parameter can be seen in Figure 10. It is interesting to note that, similarly to the cases of the 9:2 and 4:1 synodic resonant NRHOs of the Bicircular computed by Boudad et al. (2020), the Lyapunov exponents of their quasi-periodic counter parts remain approximately the same throughout the continuation procedure. Although not visible in Figure 10a, a bifurcation along one of the stable (zero) Lyapunov exponent pairs exists from the start of the perturbation, but does not depart significantly from the zero value (around 0.003 for $e = 0.0549$).

Finally, we verify that no discontinuity exists along the continuation, and as such, conclude that the resonant quasi-periodic NRHOs hereby computed represent the dynamical substitutes of the CR3BP's resonant NRHOs in the Elliptic-Circular model.

533

## 7. Torus maps for mission analysis

534

535

In the context of future missions, such as the Lunar Gateway mission, it is relevant to analyze the quasi-periodic solutions hereby presented through a mission analysis lens. To do so, we make use of two-dimensional torus maps. These maps allow us to display
<figure><img src="image_12.png" /><figcaption>Fig. 10: Lyapunov exponents plotted against the continuation parameter for each of the quasi-periodic resonant NRHOs during the continuation to the Elliptic-Circular model: (a) 9:2 synodic resonance (b) 4:1 synodic resonance and (c) 4:1 sidereal resonance.</figcaption></figure>

useful physical information over a representation of the quasi-periodic trajectories, which proves to be a significant tool for mission design and analysis.

As previously mentioned, quasi-periodic tori are characterized by their incommensurate frequencies, $\omega \in T^d$, which are constant in time. This means that by knowing the torus frequencies we can fully describe the variation of the torus angles $\theta_i$ as a function of time (Villegas-Pinto et al., 2021). For a two-dimensional torus we can write

$$ \begin{cases} \theta_0(t) = \theta_0(t_0) + \omega_0(t - t_0), \\ \theta_1(t) = \theta_1(t_0) + \omega_1(t - t_0), \end{cases} \qquad (46) $$

for some initial time $t_0$ and initial angles $\theta_i(t_0)$. This parametrization allows us to describe, based on the initial angle coordinates, how a trajectory will move on the surface of a quasi-periodic torus by means of what we call torus maps. These maps allow us to map states from $T^2$ to $R^n$ and vice-versa. This can be used to choose the best orbit insertion or departure conditions, to analyze the evolution of significant parameters such as orbit altitude, latitude, longitude, to predict eclipse events, and more. In order to obtain the torus maps, we must integrate the initial invariant circle of the quasi-periodic tori and sample the trajectories at an equally spaced grid over the two angle coordinates, $\theta_0$ and $\theta_1$. Later, we can also interpolate to approximate the values between the different nodes. Along the $\theta_0$ direction we typically choose a large number of points in order to avoid skipping sections where the dynamics are faster. In order to obtain the torus states at a regular grid we need to undo the rotation associated with the integration of the trajectory, which is done using the same rotation operator described in Eq. (16) (see Villegas-Pinto et al. (2021) for more details).

Applying this technique to the analysis of the synodic and sidereal resonant quasi-periodic NRHOs, we generate two types of torus maps. For one type, we plot the eclipses, or, more accurately, the Sun disk visibility - where a value of one corresponds to full sunlight, zero to total eclipse, and anything in between to a partial or annular eclipse - using the conic eclipse model and shadow function from Montenbruck & Gill (2000) for both Earth's and the Moon's shadows. A schematic representation of the conic eclipse model used can be seen in Figure 11, where $R_S$ represents the Sun radius and $R_i$ is the radius of either the Earth or the Moon.

For the other type of torus maps, we plot the altitude with respect to the Moon surface. These torus maps can be seen in Figures 12 to 14. On the left side of these figures we can see the Sun visibility along the torus surface, whereas the right-hand side of the figures show the altitude with respect to the Moon surface. These latter altitude torus maps are overlayed with the eclipses shown in bright green dotted lines for better visibility and analysis, albeit now without any distinction as to the type of eclipse.
<figure><img src="image_13.png" /><figcaption>Fig. 11: Schematic representation of the conic eclipse model used.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Fig. 12: Torus maps for the 9:2 synodic resonant quasi-periodic NRHO: (a) eclipse map showing fraction of Sun disk visible (b) altitude with respect to Moon surface with the eclipses overlayed in bright green dotted line.</figcaption></figure>

By analyzing these figures, we quickly verify that the eclipse avoidance properties of the synodic resonant trajectories are maintained when transitioning to Elliptic-Circular model. The maps for the 9:2 synodic resonant quasi-periodic torus seen in Figure 12 present only two small eclipse regions, close to $\theta_0 = 3\pi/4$ and $\theta_0 = 7\pi/4$, which, as expected, coincide with the $z = 0$ crossings and cover only a section of the invariant circle. In fact, it is interesting to note that without any maneuvers, it would be possible to avoid said regions for at least up to three full revolutions of $\theta_0$, i.e. $3 \times 2\pi$, by choosing appropriate initial conditions for insertion, e.g. the line starting at $\theta_0 = 0$ between $\theta_1 = \pi/2$ and $\theta_1 = 3\pi/4$ approximately. Since this quasi-periodic torus corresponds to a 9:2 synodic resonance, following such a trajectory would enable a spacecraft to avoid any type of eclipse for at least up to $3 \times 2T_{syn}$, approximately 177 days. Eventually, following the same trajectory would lead to crossings with the eclipse region, however, it would be possible to avoid these by performing correction or phasing maneuvers to transfer between two hyper-points of the torus' surface, similarly to how previous work has demonstrated transfers between intersecting tori Henry & Scheeres (2020). We also note that only Moon eclipses are detected for this resonance.

While the 9:2 synodic resonant presents two eclipse regions, we see from Figure 13 that the 4:1 synodic resonant quasi-periodic NRHO does not seem to cross any eclipse regions in the Elliptic-Circular model. This is likely due to the combination of the optimal synodic resonant geometry and the larger perilune radii visible in the 4:1 synodic resonance.
<figure><img src="image_15.png" /><figcaption>Fig. 13: Torus maps for the 4:1 synodic resonant quasi-periodic NRHO: (a) eclipse map showing fraction of Sun disk visible (b) altitude with respect to Moon surface.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Fig. 14: Torus maps for the 4:1 sidereal resonant quasi-periodic NRHO: (a) eclipse map showing fraction of Sun disk visible (b) altitude with respect to Moon surface with the eclipses overlayed in bright green dotted line.</figcaption></figure>

In contrast, the torus maps of the 4:1 sidereal quasi-periodic NRHO seen in Figure 14 display four sets of two eclipse regions, spaced approximately by $\pi/4$ along the $\theta_0$ direction and by $\pi$ along the $\theta_1$ direction. These correspond of course to the four perilune regions and to opposite sides of the Moon. Because sidereal resonances no longer have such a favorable eclipse avoidance geometry, it is normal to see more eclipse regions for this quasi-periodic torus. Nevertheless, passages through these regions could be minimized by employing the strategies described above. Additionally, as with the 9:2 synodic resonance, only Moon eclipses are detected.
Finally, it should be noted that the Elliptic-Circular model does not take into account the inclination of each of the bodies'
orbital planes. Such a consideration would be made when translating these results into a high-fidelity ephemeris model and would
likely lead to an even smaller number of eclipse events. Further analyses can be made by means of torus maps, depending on the
physical quantities we are interested in representing along the tori. These empower mission designers to choose ideal geometries
and conditions for a variety of operations, transfers, or science observations.

8. Conclusion

This paper has presented and detailed a higher-fidelity dynamical model called the Elliptic-Circular Restricted Four-Body
Problem (ECR4BP) that incorporates the eccentricity of the primary and secondary bodies' orbits and a third body in a circular
co-planar orbit about their barycenter. Motivated by the model's suitability to study the cislunar environment and by the future
exploration of the cislunar space by the Lunar Gateway and associated missions, we have investigated the quasi-periodic structures
of the Earth-Moon-Sun system that are the dynamical equivalents to periodic resonant Near-Rectilinear Halo Orbits of the Circular
Restricted Three-Body Problem. We have presented the steps and methods necessary to compute these quasi-periodic tori. We
have focused on synodic and sidereal resonant orbits, which are resonant with the period of the Earth-Moon-Sun and the period
of the Earth-Moon system, respectively. These resonant periodic orbits produce resonant two-dimensional quasi-periodic tori in
the ECR4BP model, which can be computed with the numerical continuation routines hereby outlined. We have analyzed these
solutions in terms of relevant mission parameters and their possible application to the Lunar Gateway's orbit, such as eclipses
events, altitude with respect to the Moon's surface, and the dynamical stability of these quasi-periodic trajectories. We find that the
stability of these solutions remains very close to the near-stable behavior presented by their periodic counterparts, although one of
the synodic resonances present a bifurcation along one of their stable eigenpairs during the condition procedure. We also find that
the geometry associated with synodic resonances verifies eclipse avoidance properties in the full ECR4BP model, in line with the
findings of previous works.

Acknowledgments

Daniel Villegas Pinto acknowledges the support of the Centre National d'Études Spatiales (CNES) and Airbus Defence and
Space for the funding of this research.

References

Andreu, M. A. (1998). *The quasi-bicircular problem*. Ph.D. thesis.

Assadian, N., & Pourtakdoust, S. H. (2010). On the quasi-equilibria of the BiElliptic four-body problem with non-coplanar motion of primaries. *Acta Astronautica, 66*(1-2), 45-58. URL: http://dx.doi.org/10.1016/j.actaastro.2009.05.014. doi:10.1016/j.actaastro.2009.05.014.

Baresi, N. (2017). *Spacecraft Formation Flight on Quasi-periodic Invariant Tori*. Ph.D. thesis University of Colorado Boulder.

Baresi, N., Olikara, Z. P., & Scheeres, D. J. (2018). Fully numerical methods for continuing families of quasi-periodic invariant tori in astrodynamics. *Journal of the Astronautical Sciences, 65*(2), 157-182. doi:10.1007/s40295-017-0124-6.

Betts, J. T. (2010). *Practical Methods for Optimal Control and Estimation Using Nonlinear Programming*.

Boudad, K. K., Howell, K. C., & Davis, D. C. (2020). Dynamics of synodic resonant near rectilinear halo orbits in the bicircular four-body problem. *Advances in Space Research, 66*(9), 2194-2214. URL: https://doi.org/10.1016/j.asr.2020.07.044. doi:10.1016/j.asr.2020.07.044.

Campagnola, S., Lot, M., & Newton, P. (2008). Subregions of motion and elliptic halo orbits in the elliptic restricted three-body problem. *Advances in the Astronautical Sciences, 130* PART 2, 1541-1555.

Castellà, E., & Jorba, À. (2000). On the vertical families of two-dimensional tori near the triangular points of the bicircular problem. *Celestial Mechanics and Dynamical Astronomy, 76*(1), 35-54. doi:10.1023/A:1008321605028.

Castellà, E., & Jorba, À. (2003). The lagrangian points in the real Earth-Moon system. In *International Conference on Differential Equations* (pp. 3-12). Hassel, Belgium.

Chakraborty, A., & Narayan, A. (2019). BiElliptic Restricted Four Body Problem. *Few-Body Systems, 60*(1). URL: https://doi.org/10.1007/s00601-018-1472-x. doi:10.1007/s00601-018-1472-x.

Chikazawa, T., Baresi, N., Campagnola, S., Ozaki, N., & Kawakatsu, Y. (2021). Minimizing eclipses via synodic resonant orbits with applications to EQUULEUS and MMX. *Acta Astronautica*, (180), 679-692. URL: https://doi.org/10.1016/j.actaastro.2020.12.028. doi:10.1016/j.actaastro.2020.12.028.
Davis, D. C., Bhatt, S. A., Howell, K. C., Jang, J. W., Whitley, R. J., Clark, F. D., Guzzetti, D., Zimovan, E. M., & Barton, G. H. (2017). Orbit maintenance and navigation of human spacecraft at cislunar near rectilinear halo orbits. *Advances in the Astronautical Sciences*, **160**, 2257–2276.

Gómez, G., & Mondelo, J. M. (2001). The dynamics around the collinear equilibrium points of the RTBP. *Physica D: Nonlinear Phenomena*, **157**(4), 283–321.
doi:10.1016/S0167-2789(01)00312-8.

Gómez, G., Simó, C., Llibre, J., & Martínez, R. (2001). *Dynamics and mission design near libration points - Vol. II fundamentals: the case of the triangular libration points* volume 3. World Scientific.

Henderson, M. E. (2002). Multiple parameter continuation: Computing implicitly defined k-manifolds. *International Journal of Bifurcation and Chaos in Applied Sciences and Engineering*, **12**(3), 451–476. doi:10.1142/S0218127402004498.

Henry, D. B., & Scheeres, D. J. (2020). Transfers between intersecting quasi-periodic tori. In *AIAA SciTech Forum* (pp. 1–17). Orlando, Florida.

Howell, K. C., & Breakwell, J. V. (1984). Almost rectilinear halo orbits. *Celestial Mechanics*, **32**(1), 29–52. doi:10.1007/BF01358402.

IAU Division Working Group (2022). Numerical Standards for Fundamental Astronomy - Astronomical Constants: Current Best Estimates (CBEs). URL: https://iau-a3.gitlab.io/NSFA/NSFA[_]cbe.html.

Jenson, E. L., & Scheeres, D. J. (2021). Semianalytical Measures of Nonlinearity Based on Tensor Eigenpairs. *AAS/AIAA Astrodynamics Specialist Conference*, (pp. 1–17).

Jorba, À., & Villanueva, J. (1997). On the persistence of lower dimensional invariant tori under quasi-periodic perturbations. *Journal of Nonlinear Science*, **7**(5), 427–473. doi:10.1007/s00329900036.

Jorba-Cuscó, M., Farrés, A., & Jorba, À. (2018). Two periodic models for the Earth-Moon system. *Frontiers in Applied Mathematics and Statistics*, 4(July), 1–14. doi:10.3389/fams.2018.00032.

Mamajek, E. E., Prsa, A., Torres, G., Harmanec, P., Asplund, M., Bennett, P. D., Capitaine, N., Christensen-Dalsgaard, J., Depagne, E., Folkner, W. M., Haberreiter, M., Hekker, S., Hilton, J. L., Kostov, V., Kurtz, D. W., Laskar, J., Mason, B. D., Milone, E. F., Montgomery, M. M., Richards, M. T., Schou, J., & Stewart, S. G. (2015). IAU 2015 Resolution B3 on Recommended Nominal Conversion Constants for Selected Solar and Planetary Properties, (pp. 1–6). URL: http://arxiv.org/abs/1510.07674. arXiv:1510.07674.

McCarthy, B., & Howell, K. C. (2022). Characterization of Families of Low-Energy Transfers to Cislunar Four-Body Quasi-Periodic Orbits. In *AIAA SciTech Forum* (pp. 1–17). San Diego, California. doi:10.2514/6.2022-1889.

Montenbruck, O., & Gill, E. (2000). *Satellite Orbits - Models, Methods, and Applications*. Springer, Berlin, Heidelberg. doi:https://doi.org/10.1007/978-3-642-58351-3.

Muñoz-Almaraz, F. J., Freire, E., Galán, J., Doedel, E., & Vanderbauwhede, A. (2003). Continuation of periodic orbits in conservative and Hamiltonian systems. *Physica D: Nonlinear Phenomena*, **181**(1-2), 1–38. doi:10.1016/S0167-2789(03)00097-6.

NASA (2019). White Paper: Gateway Destination Orbit Model: A Continuous 15 Year NRHO Reference Trajectory, URL: https://ntrs.nasa.gov/citations/20190030294.

NASA (2020). NASA's *Lunar Exploration Program Overview*. Technical Report September https://www.nasa.gov/sites/default/files/atoms/files/artemis_plan-20200921.pdf. URL: https://www.nasa.gov/sites/default/files/atoms/files/artemis[_]plan-20200921.pdf.

Olikara, Z. P. (2016). *Computation of quasi-periodic tori and heteroclinic connections in astrodynamics using collocation techniques*. Ph.D. thesis University of Colorado Boulder.

Olikara, Z. P., & Scheeres, D. J. (2012). Numerical method for computing quasi-periodic orbits and their stability in the restricted three-body problem. *Advances in the Astronautical Sciences*, **145**, 911–930.

Oshima, K. (2022). Continuation and stationkeeping analyses on planar retrograde periodic orbits around the Earth. *Advances in Space Research*, **69**(5), 2210–2222. URL: https://doi.org/10.1016/j.asr.2021.12.020. doi:10.1016/j.asr.2021.12.020.

Rosales, J. J., Jorba, À., & Jorba-Cuscó, M. (2021). Families of Halo-like invariant tori around L2 in the Earth-Moon Bicircular Problem. *Celestial Mechanics and Dynamical Astronomy*, **133**(4), 1–30. URL: https://doi.org/10.1007/s10569-021-10012-0. doi:10.1007/s10569-021-10012-0.

Scheeres, D. J. (1998). The Restricted Hill Four-Body Problem with Applications to the Earth-Moon-Sun System. *Celestial Mechanics and Dynamical Astronomy*, **70**(2), 75–98. doi:10.1023/A:1026498608950.

Scheeres, D. J. (2012). *Orbital Motion in Strongly Perturbed Environments*. doi:10.1007/978-3-642-03256-1.

Scheeres, D. J., & Bellerose, J. (2005). The Restricted Hill Full 4-Body problem: Application to spacecraft motion about binary asteroids. *Dynamical Systems*, **20**(1), 23–44. doi:10.1080/1468936042000281321.

Scheuerle, S. T., & Howell, K. C. (2021). Characteristics and Analysis of Families of Low-Energy Ballistic Lunar Transfers. In *AAS/AIAA Astrodynamics Specialist Conference* (pp. 1–20). Big Sky, Montana (Virtual).

Seydel, R. (2009). *Practical Bifurcation and Stability Analysis* volume 5. (3rd ed.). Springer Science & Business Media. URL: http://www.springerlink.com/index/10.1007/978-0-387-75847-3. doi:10.1007/978-0-387-75847-3.

Szebehely, V. (1967). *Theory of orbits*. Academic Press Inc. doi:10.1016/0016-0032(67)90217-7.

Villegas-Pinto, D., Baresi, N., Hestroffer, D., & Canalias, E. (2021). On the numerical computation of quasi-periodic families and applications to the Martian Moons Exploration mission. In *ICATT 2021*.

Wieczorek, M. A., Jolliff, B. L., Khan, A., Pritchard, M. E., Weiss, B. P., Williams, J. G., Hood, L. L., Righter, K., Neal, C. R., Shearer, C. K., McCallum, I. S., Tompkins, S., Hawke, B. R., Peterson, C., Gillis, J. J., & Bussey, B. (2006). The constitution and structure of the Lunar interior. *Reviews in Mineralogy and Geochemistry*, **60**, 221–364. doi:10.2138/rmg.2006.60.3.

Zimovan, E. M., Howell, K. C., & Davis, D. C. (2017). Near rectilinear halo orbits and their application in cis-lunar space. In *3rd IAA Conference on Dynamics and Control of Space Systems*. Moscow, Russia.

Wieczorek, M. A., Jolliff, B. L., Khan, A., Pritchard, M. E., Weiss, B. P., Williams, J. G., Hood, L. L., Righter, K., Neal, C. R., Shearer, C. K., McCallum, I. S., Tompkins, S., Hawke, B. R., Peterson, C., Gillis, J. J., & Bussey, B. (2006). The constitution and structure of the Lunar interior. *Reviews in Mineralogy and Geochemistry*, **60**, 221–364. doi:10.2138/rmg.2006.60.3.

Zimovan, E. M., Howell, K. C., & Davis, D. C. (2017). Near rectilinear halo orbits and their application in cis-lunar space. In *3rd IAA Conference on Dynamics and Control of Space Systems*. Moscow, Russia.