Interplanetary transfers using stable and unstable
manifold tubes originating at L1 and L2 in the elliptic
restricted three-body problems

Erica Scantamburlo, Massimiliano Guzzo, Rocio Isabel Paez

► To cite this version:

Erica Scantamburlo, Massimiliano Guzzo, Rocio Isabel Paez. Interplanetary transfers using stable and unstable manifold tubes originating at L1 and L2 in the elliptic restricted three-body problems. 2021. hal-03436121

HAL Id: hal-03436121

https://hal.science/hal-03436121

Preprint submitted on 19 Nov 2021

**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
# Interplanetary transfers using stable and unstable manifold tubes originating at $L_1$ and $L_2$ in the elliptic restricted three-body problems

Erica Scantamburlo¹, Massimiliano Guzzo², and Rocio Isabel Paez³

¹Centro di Ateneo di Studi e Attività Spaziali “Giuseppe Colombo” - CISAS, Università degli Studi di Padova

²,³Dipartimento di Matematica “Tullio Levi Civita”, Università degli Studi di Padova

November 19, 2021

## Abstract

The solutions originating at the Lagrangian points $L_1$ and $L_2$ of the restricted three-body problem are used in space flight dynamics to find useful station-keeping orbits as well as convenient low-energy transfers between planets and satellites of the Solar System. The circular restricted three-body problem (CR3BP) provides the simplest model where the Lagrangian points, as well as the related stable and unstable manifold tubes, are defined. Nevertheless, their use for Solar System applications requires to consider hierarchical extensions of the model. The elliptic restricted three-body problem (ER3BP) is the most important extension to consider for space flight dynamics. In this paper we discuss the differences of the manifold tubes of the ER3BP with respect to the CR3BP which have an impact in the transfers between orbits in the manifold tubes related to different planets. We find that there is a threshold value for the eccentricity of the planet, depending on its reduced mass, which changes drastically the distribution of the longitudes of the first perycenter (or apocenter) of the orbits in the manifold tubes leaving or approaching the Lagrangian points $L_1, L_2$. For example, such a difference changes the design of the Hohmann transfers between two orbits of the manifold tubes when the planets Mercury or Mars are concerned.

## 1 Introduction

In the last decades, the solutions originating at the Lagrangian points $L_1$ and $L_2$ of the circular restricted three-body problem (CR3BP) have been widely analyzed in connection to transport phenomena, close encounters of comets and asteroids as well as low-energy transfer orbits (for example, see Jorba and Masdemont, 1999; Palacián and Yanguas, 2004; Ren et al., 2012; Guzzo and Lega, 2013, 2015, 2017, 2018; Paez and Guzzo, 2020; Ceccaroni et al., 2016; Pucacco, 2019 and Simó, 1999; Gómez et al., 2000, 2004; Koon et al., 2008; Zanzottera et al., 2012; Pergola et al., 2009; Belbruno et al., 2010; Topputo et al., 2005; Parker and Anderson, 2013; Farquhar et al., 2004; Howell et al., 2006; Howell and Kakoi, 2006). As a matter of fact several space missions (such as ISEE-3, Wind, SOHO, WMAP, Genesis, Herschel and Planck) were designed by exploiting the Sun-Earth Lagrangian points $L_1$ and $L_2$, and other planned future missions will reach the vicinity of such Lagrangian points. If we consider the Sun-Mercury or the Sun-Mars systems, the CR3BP provides only a crude approximation of the dynamics; as a matter of fact the eccentricities of the two planets are currently about 0.2 and 0.09; for this reason, if we plan to apply the space manifold dynamics
between the Sun-Earth to the Sun-Mercury or Sun-Mars Lagrangian points, a deep analysis of the solutions originating at $L_1$ and $L_2$ in the elliptic restricted three-body problem (ER3BP) is required. When we consider a generalization of space manifold dynamics to the ER3BP several difficulties appear which on the one hand are of substantial nature (such as the lack of a constant of motion generalizing the Jacobi integral), on the other hand are of technical nature (the development of perturbation theories is complicated by the implementation of a partially hyperbolic Floquet theory). In Paez and Guzzo (2021) a generalization of the space manifold dynamics to the solutions originating at $L_1, L_2$ of the ER3BP was done through the construction of Hamiltonian Floquet-Birkhoff normal forms.

The aim of the paper is to analyze the differences of the manifold tubes originating at the Lagrangian solutions in the planar ER3BP and the CR3BP, and to discuss the impact they have in the design of interplanetary transfers. Using the definition and the computation of the manifold tubes provided in Paez and Guzzo (2021), we first characterize the evolution of the orbital elements of a spacecraft whose motion is described by a solution of the ER3BP that belong to the manifold tubes of a generic Sun-$P_j$ system. Specifically, we consider any solution whose initial condition belongs to the local unstable manifold tube originating at $L_1$ or $L_2$. For an initial time-span the solution remains close to the Lagrangian one, oscillating around it. Then, the spacecraft exits from the planet's Hill sphere and we follow its evolution until it arrives sufficiently far from the planet so that its dynamics is well approximated by the Sun-spacecraft two-body problem. After that moment, the solution is well approximated by a Keplerian one up to an eventual additional close encounter with $P_j$, and therefore the semi-major axis, eccentricity and longitude of perihelion of the spacecraft remain approximately constant. The solutions whose initial condition belong to the local stable manifold tube are characterized by the same behavior with reference to negative time-spans. We denote by $K_{j,L_i}^s$ and $K_{j,L_i}^u$ the family of all these Keplerian orbits which are obtained for all the possible initial conditions in the local stable or unstable manifold tubes originating at $L_i$ of the Sun-$P_j$ ER3BP, and for all the values of initial true anomaly of the planet. The idea is to exploit the orbits in the sets $K_{h,L_i}^{s/u}$ and $K_{k,L_m}^{s/u}$ for the design of an interplanetary transfer between the two planets $P_h$ and $P_k$.

With evidence, the construction of the families $K_{j,L_i}^{s/u}$ can be defined using the simpler CR3BP, and one may argue if indeed the more complicated ER3BP is worthy to be considered. Here we find the main differences between the two models: depending on the value of the eccentricity $e_j$ of the orbit of planet $P_j$, we find a big difference in the distribution of the orbital parameters of the Keplerian orbits in the sets $K_{j,L_i}^{s/u}$ when computed using the ER3BP rather than the CR3BP. Precisely, we find that for:

$$e_j \gtrsim \frac{3r_{\text{Hill}}}{a_j}, \qquad (1)$$

where $a_j, e_j$ denote the semi-major axis and eccentricity of the planet $P_j$ and $r_{\text{Hill}}$ the Hill radius, the distribution of the longitudes of perihelion of the orbits in the above families is bounded; conversely, for $e_j \lesssim (3r_{\text{Hill}}/a_j)$ the distribution is spread in $[0, 2\pi]$ (the latter case, for $e_j = 0$, corresponds to the CR3BP). From a space flight dynamics point of view, this implies that:

- if both the eccentricities of the planets $P_h$ and $P_k$ are large (i.e. inequality (1) is satisfied for $P_h, P_k$), then the distribution of longitudes of the perihelion for the orbits in the families $K_{h,L_i}^{s/u}$ and $K_{k,L_m}^{s/u}$ of each planet is bounded in some intervals, and if such intervals do not match it is not possible to perform an Hohmann transfer<fnref n="1" /> between $P_h$ and $P_k$ via the families

<footnote n="1">¹Since we find that all the orbits in the families have eccentricity different from 0 (see Sections 2 and 3 for details, and also Fig. 7), here we need to refer to the extended definition of Hohmann transfers between co-axial elliptic orbits (see, for example, Curtis, 2014). Precisely a transfer between two different co-axial elliptic orbits is obtained with two</footnote>
$K_{h,L_i}^{s/u}$, $K_{k,L_m}^{s/u}$; in such a case, other types of transfers may be considered, in the context of the ER3BP.

- if both the eccentricities of the planets $P_h$ and $P_k$ are small (i.e. inequality (1) is not satisfied for $P_h, P_k$), then the distribution of longitudes of the perihelion for the orbits in the families $K_{h,L_i}^s, K_{k,L_m}^s$ of each planet is distributed in $[0, 2\pi]$, and we have an abundance of co-axial orbits in the families which are suitable for an Hohmann transfer. This is the situation where the CR3BP is a suitable model for the transfer.

- the most interesting situation arises when one planet satisfies the inequality (1), while the other does not. In such a case the longitudes of perihelion of the orbits in the families $K_{h,L_i}^{s/u}$ of the latter planet spread in $[0, 2\pi]$, and therefore we still find an abundance of co-axial orbits suitable for the Hohmann transfer. Nevertheless, since the longitudes of perihelion of the orbits in the families $K_{k,L_m}^{s/u}$ of the former planet are in a bounded interval, it is necessary to properly synchronize the maneuvers of the transfer in suitable time windows. This is the situation where the ER3BP must be considered for the transfer.

We apply these results to the Sun–Mercury, Sun–Earth and Sun–Mars systems: since both Mercury and Mars satisfy inequality (1), while the Earth (as well as Venus) does not fulfill such condition, the distribution of the longitudes of perihelion is bounded in the Mercury and Mars systems, and is spread in $[0, 2\pi]$ in the Earth system. This implies that it is possible to determine an Hohmann transfer between the manifold tubes originating at a Lagrangian points of the Sun–Earth and Sun–Mercury, or Sun–Mars ER3BP. Then, we calculate the transfer between the Sun–Earth $L_1$ and the Sun–Mercury $L_2$, and between the Sun–Earth $L_2$ and the Sun–Mars $L_1$ points using the decoupled ER3BPs<fnref n="2" />. Finally we check that the transfers, whose initial conditions are analytically computed using two decoupled ER3BPs, remain effective also in the bi-elliptic four body problem (BE4BP), namely a model in which a spacecraft is gravitationally attracted by the Sun and the two planets that are performing an elliptic motion around the Sun, for the whole considered time interval.

The paper is organized as follow: in Sec. 2 we first define the planar BER4BP, and then we briefly summarize the procedure used in Paez and Guzzo (2021) for the computation of the solutions originating at the Lagrangian points $L_1$ and $L_2$ in the ER3BP. Then we analyze the orbital elements evolution of the manifold tubes. Finally, in Sec. 3 we apply the results found to an interplanetary transfer between the Earth and Mercury, and the Earth and Mars systems using both the decoupled ER3BPs and the BER4BP.

# 2 Study of the spacecraft orbital elements around $L_1$ and $L_2$ in the elliptic restricted three-body problem

## 2.1 The planar bi-elliptic restricted four-body problem

To model the transfer of a spacecraft between the manifold tubes of the planets $P_h, P_k$ we first simplify the equations of motion in the full Solar System to the planar bi-elliptic restricted four-body problem (BER4BP). Denoting by $\mathbf{r}_j$ and $\mathbf{r}$ the position vectors of the $j$-th planet and of the spacecraft in a heliocentric reference frame $Oxyz$, by $\mu$ and $\epsilon_j$ the gravitational parameters of the

<footnote>impulsive maneuvers: the first maneuver is obtained by applying a tangential $\Delta v$ at a perycenter/apocenter of the first orbit; the spacecraft is thus injected into the transfer orbit and when it arrives to the next perycenter/apocenter it is injected into the final keplerian orbit with an additional tangential $\Delta v$.</footnote>

<footnote n="2">²We do not compute transfers between the Sun–Earth and Sun–Venus Lagrangian points since both planets do not satisfy inequality (1), and therefore the transfers are suitably described with the CR3BP.</footnote>
Sun and of the j-th planet, and by $(a_j, e_j, \varpi_j, f_j)$ the orbital elements of planet $P_j$, the equations of motion of the spacecraft are

$$ \ddot{\mathbf{r}} = -\mu \frac{\mathbf{r}}{\|\mathbf{r}\|^3} - \sum_{j=h,k} \varepsilon_j \left( \frac{\mathbf{r}_j}{\|\mathbf{r}_j\|^3} + \frac{\mathbf{r} - \mathbf{r}_j}{\|\mathbf{r} - \mathbf{r}_j\|^3} \right). \quad (2) $$

Then, we compare the dynamics of the BER4BP with that of the two decoupled Sun-$P_h$ and Sun-$P_k$ elliptic restricted three-body problems (ER3BPs). For this reason, for each planet we introduce the Sun-planet heliocentric rotating-pulsating reference frame OXY which is centered at the Sun, and whose axes are parallel to the unit vectors

$$ \mathbf{u}_{r,j} = \frac{\mathbf{r}_j}{\|\mathbf{r}_j\|}, \quad \mathbf{u}_{t,j} = \frac{\mathbf{r}_j \times \mathbf{u}_{n,j}}{\|\mathbf{r}_j \times \mathbf{u}_{n,j}\|}, \quad (3) $$

where $\mathbf{u}_{n,j} = \frac{\dot{\mathbf{r}}_j \times \mathbf{r}_j}{\|\dot{\mathbf{r}}_j \times \mathbf{r}_j\|}$ and $j = h, k$; in such a reference frame, the planet $P_j$ is located on the X axis.

We investigate the dynamics close to the Lagrangian solutions of the Sun-$P_j$ ER3BP using the Floquet-Birkhoff normal forms defined in Paez and Guzzo (2021).

## 2.2 The Floquet-Birkhoff normal form in the elliptic restricted three-body problem

It is usual to describe the planar ER3BP using a rotating-pulsating reference frame and standard Sun-planet units of measure (so that the sum of Sun and planet gravitational parameters is one, the semi-major axis of the planet orbit is one, and the revolution period of the planet is $2\pi$). The Hamilton function of the ER3BP in the heliocentric rotating-pulsating reference frame is:

$$ \mathcal{H}(X, Y, P_X, P_Y, f_j; \varepsilon_j, e_j) = \frac{1}{2} \left[ (P_X + Y)^2 + (P_Y - X)^2 \right] - \frac{1}{1 + e_j \cos(f_j)} \left[ \frac{1}{2}(X^2 + Y^2) + \frac{1}{\sqrt{X^2 + Y^2}} \right] \\ + \frac{\varepsilon_j}{1 + e_j \cos(f_j)} \left( X + \frac{1}{\sqrt{X^2 + Y^2}} - \frac{1}{\sqrt{(X-1)^2 + Y^2}} \right) \quad (4) $$

where $\varepsilon_j$ denotes the gravitational parameter of the planet, and $P_X$ and $P_Y$ represent the conjugate momenta to X and Y respectively. We remark that the Hamiltonian is not autonomous since it depends explicitly on the planet true anomaly $f_j$, which is also the independent variable of Hamilton's equations:

$$ \frac{dX}{df_j} = \frac{\partial \mathcal{H}}{\partial P_X}, \quad \frac{dY}{df_j} = \frac{\partial \mathcal{H}}{\partial P_Y}, \quad \frac{dP_X}{df_j} = -\frac{\partial \mathcal{H}}{\partial X}, \quad \frac{dP_Y}{df_j} = -\frac{\partial \mathcal{H}}{\partial Y}. \quad (5) $$

The advantage of the rotating-pulsating reference frame formulation if that the Hamiltonian system (5) has five equilibrium points $L_1, \dots, L_5$ in the same location of the Lagrangian points computed in the Sun-planet CR3BP. We focus our attention on the two collinear equilibria $L_1$ and $L_2$ (located at $(X, Y, P_X, P_Y) = (X_{L_1}, 0, 0, X_{L_2})$, $i = 1, 2$). Following Paez and Guzzo (2021) we construct the manifold tubes originating at $L_1, L_2$ of the ER3BP using a suitable Floquet-Birkhoff normal form. We here summarize the construction of such normal forms.

Let us first rewrite the Hamiltonian (4) using the translated Cartesian variables:

$$ q_1 := X - X_{L_1}, \quad q_2 := Y, \quad p_1 := P_X, \quad p_2 := P_Y - X_{L_1} \quad (6) $$
and then we expand it in Taylor series at the origin:

$$H(\mathbf{q}, \mathbf{p}; f_j) := \sum_{\alpha=2}^{\infty} H_{\alpha}(\mathbf{q}, \mathbf{p}; f_j) \qquad (7)$$

where $H_\alpha(\mathbf{q}, \mathbf{p}; f_j)$ is a polynomial of degree $\alpha$ in $(\mathbf{q}, \mathbf{p})$. In the paper Paez and Guzzo (2021), a canonical transformation:

$$(\mathbf{q}, \mathbf{p}; f_j) = \mathcal{C}(\mathbf{Q}, \mathbf{P}; f_j) \qquad (8)$$

has been constructed conjugating the Hamiltonian (7) to:

$$K(\mathbf{Q}, \mathbf{P}; f_j) = \sum_{\nu=1}^{N/2} K_{2\nu}(\mathbf{Q}, \mathbf{P}) + R_{N+1}(\mathbf{Q}, \mathbf{P}; f_j) \qquad (9)$$

where $K_{2\nu}(\mathbf{Q}, \mathbf{P})$ denotes an autonomous polynomial of degree $2\nu$ depending on $(\mathbf{Q}, \mathbf{P})$ only through the combinations of the terms $\frac{Q_1^2+P_1^2}{2}$ and $Q_2P_2$; the remainder $R_{N+1}(\mathbf{Q}, \mathbf{P}; f_j)$ is a polynomial of order $N+1$, possibly dependent on $f_j$; $N$ is such that the frequencies of the problem do not satisfy resonances of order smaller than $N$ (see the paper Paez and Guzzo, 2021, for details). Since the transformation $\mathcal{C}(\mathbf{Q}, \mathbf{P}; f_j)$ is obtained from the composition of a Floquet transformation and $N-2$ Birkhoff transformations, the Hamiltonian (9) has been called Floquet-Birkhoff normal form.

A classification of the orbits in a neighborhood of the equilibrium $L_i$ is provided by the flow of the integrable Hamiltonian:

$$\mathcal{K}(\mathbf{Q}, \mathbf{P}) := \sum_{\nu=1}^{N/2} K_{2\nu}(\mathbf{Q}, \mathbf{P}), \qquad (10)$$

which is obtained from (9) by neglecting the remainder $R_{N+1}(\mathbf{Q}, \mathbf{P}; f_j)$. We emphasize that when the variables $(Q_1, P_1)$ represent the elliptic motion, while $(Q_2, P_2)$ the hyperbolic one, the lowest order approximation of $\mathcal{K}$ is

$$K_2(\mathbf{Q}, \mathbf{P}) = \sigma \frac{Q_1^2 + P_1^2}{2} + \lambda Q_2 P_2 \quad \sigma, \lambda > 0. \qquad (11)$$

The flow of the approximated normal form (11) has a family of periodic orbits, the planar Lyapunov orbits identified by $Q_2, P_2 = 0$ and $(Q_1, P_1) = \sqrt{2I_1}(\sin \phi_1, \cos \phi_1)$, with $\phi_1 \in [0, 2\pi]$ and $I_1 > 0$ suitably small. These Lyapunov orbits have local stable manifold identified by $(Q_1, P_1) = \sqrt{2I_1}(\sin \phi_1, \cos \phi_1)$, $P_2 \neq 0$ suitably small and $Q_2 = 0$, while the local unstable manifold is identified by $(Q_1, P_1) = \sqrt{2I_1}(\sin \phi_1, \cos \phi_1)$, $Q_2 \neq 0$ suitably small and $P_2 = 0$.

When we map back these sets to the Cartesian rotating-pulsating reference frame with the canonical transformation

$$(X, Y, P_X, P_Y; f_j) = (X_{L_i}, 0, 0, X_{L_i}; 0) + \mathcal{C}(\mathbf{Q}, \mathbf{P}; f_j)$$

we obtain:

- the planar Lyapunov tori:

$$\mathcal{T}_{I_1} = \bigcup_{f_j, \phi_1 \in [0, 2\pi]} \{ (X, Y, P_X, P_Y; f_j) = (X_{L_i}, 0, 0, X_{L_i}; 0) + \mathcal{C}(\mathbf{Q}, \mathbf{P}; f_j) : \\ (Q_1, P_1) = \sqrt{2I_1}(\sin \phi_1, \cos \phi_1), Q_2, P_2 = 0 \}; \qquad (12)$$
- the local unstable tubes of $\mathcal{T}_{I_1}$:

$$ \begin{aligned} W_{I_1}^{u_{loc}} &= \bigcup_{f_j, \phi_1 \in [0, 2\pi]} \{ (X, Y, P_X, P_Y; f_j) = (X_{L_i}, 0, 0, X_{L_i}; 0) + C(Q, P; f_j) : \\ & \qquad (Q_1, P_1) = \sqrt{2I_1}(\sin \phi_1, \cos \phi_1), |Q_2| \in (0, \epsilon), P_2 = 0 \}, \end{aligned} \quad (13) $$

where $\epsilon$ is such that $(Q, P)$ is in the domain of the transformation $C(Q, P; f_j)$;

- the local stable tubes of $\mathcal{T}_{I_1}$:

$$ \begin{aligned} W_{I_1}^{s_{loc}} &= \bigcup_{f_j, \phi_1 \in [0, 2\pi]} \{ (X, Y, P_X, P_Y; f_j) = (X_{L_i}, 0, 0, X_{L_i}; 0) + C(Q, P; f_j) : \\ & \qquad (Q_1, P_1) = \sqrt{2I_1}(\sin \phi_1, \cos \phi_1), |P_2| \in (0, \epsilon), Q_2 = 0 \}, \end{aligned} \quad (14) $$

where $\epsilon$ is such that $(Q, P)$ is in the domain of the transformation $C(Q, P; f_j)$.

Then, we define the global stable and unstable manifolds $W_{I_1}^s, W_{I_1}^u$ of the planar Lyapunov tori by propagating the initial conditions in the sets $W_{I_1}^{s_{loc}}, W_{I_1}^{u_{loc}}$ with the flow of the ER3BP.

**Remark.** We remark that the Lyapunov torus is indeed invariant, and the solutions in $W_{I_1}^s, W_{I_1}^u$ converge (in the future or in the past) to it only in the approximation which is obtained by neglecting the remainder of the Floquet-Birkhoff normal form $R_{N+1}$. The errors introduced by such an approximation are negligible for space flight dynamics, as soon $N$ is suitably large (we will use $N = 6$).

In Fig. 1 we show the projection of the manifold tubes originating at $L_1$ and $L_2$ in the Sun–Earth and Sun–Mercury ER3BP (first four panels) on the plane $OXY$; the initial conditions have been chosen in the local stable and unstable manifolds $W_{loc}^s, W_{loc}^u$ with $f_{j0} = 0$, and then are propagated with the flow of the ER3BP till the orbits arrive at the first aphelion (or perihelion respectively) which is sufficiently far from the planet. We note that the manifold tube in the Sun–Mercury system enlarges when the solution $(X(f_j), Y(f_j))$ is getting away from the planet (panel d)); this behavior is less evident for the Sun–Earth system (panel c)). This is an effect introduced by the rotating-pulsating reference frame; as a matter of fact the representation of the same manifold tube in the rescaled variables $\frac{1-e_j^2}{1+e_j \cos f}(X(f_j), Y(f_j))$ (panel f)) are not affected by this enlargement, while we appreciate the oscillation due to the motion of the Lagrangian solutions in the non-pulsating reference frame.

## 2.3 Orbital element representation of the manifold tubes of the ER3BP in the inertial reference frame

In this subsection we provide evidence that there exists a threshold value for the eccentricity and the gravitational parameter of the planet $P_j$ which changes drastically the distribution of the longitude of perihelion of the orbits in the manifold tubes originating at $L_1$ and $L_2$ in the ER3BP. For example, if the Sun–$P_h$ and the Sun–$P_k$ systems fulfill such threshold value, and the distributions of the longitude of perihelion do not match in the two ER3BPs, this implies that a Hohmann transfer can not be used to transfer a spacecraft from the unstable manifold tube originating at a collinear Lagrangian point of $P_h$ to the stable manifold tube originating at a Lagrangian point of $P_k$. Conversely, Hohmann transfer exists when at least one of the planets has eccentricity below the threshold value.

For this reason we study the representation of the manifold tubes of both Sun–$P_h$ and Sun–$P_k$ ER3BPs using the orbital elements computed in the same heliocentric (non-rotating and non-pulsating) reference frame. The transformation from the rotating-pulsating reference frame of
<figure><img src="image_2.png" /><figcaption>Figure 1: Panels a) and c) show the unstable manifold tube originating at the Sun–Earth $L_1$ for $I_1 = 5 \times 10^{-7}$, $Q_2 = 10^{-6}$, $f_{j0} = 0$ and $\phi_1 \in [0, 2\pi]$, while panels b) and d) show the stable manifold tube originating at the Sun–Mercury $L_2$ for $I_1 = 5 \times 10^{-7}$, $P_2 = 10^{-6}$, $f_{j0} = 0$ and $\phi_1 \in [0, 2\pi]$. The black dots belong to the orbits in the manifold tubes up to the stopping points (aphelion for the Sun–Earth system, perihelion for the Sun–Mercury system), while the gray dots belong to the same orbits, but for a later time. The red points indicate the location of the Lagrangian point. Panels e) and f) display the same manifolds using the rescaled coordinates associated to the Sun–Earth and Sun–Mercury systems respectively (the reference frame is still a rotating, but not pulsating, one).</figcaption></figure>
planet $P_j$ to such a heliocentric (non-rotating and non-pulsating) reference frame is denoted by:

$$ (x, y, v_x, v_y) = T_j(X, Y, P_X, P_Y; f_j) := T(X, Y, P_X, P_Y; a_j, e_j, f_j, \varpi_j) \quad (15) $$

where $X, Y, P_X, P_Y$ represent the canonical coordinates and momenta of the rotating-pulsating reference frame; $v_x := dx/dt$, $v_y := dy/dt$ denote the Cartesian components of the velocity in the heliocentric reference frame; the dependence on the planet $P_j$ is absorbed by $a_j, e_j, f_j, \varpi_j$.

**A family of orbits in the manifold tubes.** Let us consider a solution $(X(f_j), Y(f_j), P_X(f_j), P_Y(f_j), f_j)$ with initial condition in $W_{I_1}^u$, with $(X(f_{j0}), Y(f_{j0}), P_X(f_{j0}), P_Y(f_{j0}), f_{j0})$ close to the torus $T_{I_1}$ obtained using the Floquet-Birkhoff normal form. We highlight the evolution of the solution in different time intervals:

(i) For an initial interval $\Delta f_1 = [f_0, f_1]$ the solution remains close to the Lagrangian point, and oscillates around it. In this time interval the representation of the solution on the *Oxy* reference frame using the transformation $T_j$ (see (15)) is characterized by orbital elements oscillating close to the orbital elements of the Lagrangian solution.

(ii) In a following interval $\Delta f_2 = [f_1, f_2]$ the solution exits from the Hill sphere of the planet abandoning the local unstable manifold, and arrives to a distance from the planet where the spacecraft-Sun two body problem provides a good approximation of the dynamics. The stopping point $f_2$ is easily identified since the time evolution of the orbital elements of the spacecraft become almost constant before $f_2$. Moreover, we choose $f_2$ to be a perihelion or aphelion (depending on the transfer).

(iii) For a long time interval $\Delta f_3 = [f_2, f_3]$, where the solution does not approach again the Hill sphere of the planet, the orbit of the spacecraft is well approximated by a Keplerian one, with orbital elements determined at $f_2$. The family $K_{I_1}^u$ of all these Keplerian orbits, obtained for all the possible initial conditions in the local unstable manifold, can be used for the design of transfers.

The same arguments apply with minor changes to the solutions in the stable manifolds $W_{I_1}^s$. We therefore provide a representation of the orbital elements for the Keplerian orbits in the families $K_{I_1}^u, K_{I_1}^s$. Depending on the values of the action $I_1$ and of the phase $\phi_1$ which determine the initial conditions on the local manifolds, and on the value of the true anomaly of the planet $f_j$, we obtain the distribution of the orbital elements of the spacecraft $a, e, f, \varpi$. We find that for small values of the eccentricity (and therefore including the CR3BP approximation), the longitudes of perihelia in the family spans the full interval $[0, 2\pi]$, thus allowing a simple synchronization of orbital transfers between the manifolds of two planets and a wide choice of orbital maneuvers. The distribution of $\varpi$ in the family becomes a smaller interval for larger values of the eccentricity, thus restricting the possibility of orbital maneuvers and, as a minor consequence, the optimizations of the $\Delta v$. The Earth and Venus fulfill the first case, while Mars and Mercury fulfill the second one.

**Orbital elements evolution of the Lagrangian solutions.** Let us first compute a convenient approximation for the evolution of the orbital elements $a, e, f, \varpi$ of the Lagrangian solutions $L_1$ and $L_2$ in the *Oxy* reference frame, whose orbits is an ellipse described by $\mathbf{r}(t) = X_{L_1}\mathbf{r}_j(t)$, where $\mathbf{r}$ and $\mathbf{r}_j$ denote the position vector of the Lagrangian solution and of the planet $P_j$ respectively in the *Oxy* reference frame. By solving $\mathbf{r}(t) = X_{L_1}\mathbf{r}_j(t)$ for all $t$, we find immediately

$$ f + \varpi = f_j + \varpi_j. \quad (16) $$
Moreover, by neglecting the difference in the gravitational parameter of $P_j$ with respect to the gravitational parameter of the Sun, we obtain:

$$a \approx \frac{a_j(e_j^2 - 1)X_{L_i}}{2e_j(X_{L_i}^3 - 1)\cos f_j + (1 + e_j^2)X_{L_i}^3 - 2} \quad (17)$$

and, by defining

$$\xi := e \cos f, \quad \eta := e \sin f, \quad (18)$$

we get

$$\xi = X_{L_i}^3 - 1 + e_j X_{L_i}^3 \cos f_j, \quad \eta = e_j X_{L_i}^3 \sin f_j; \quad (19)$$

i.e. $(\xi, \eta)$ describes a circumference parameterized by the planet true anomaly $f_j$, whose radius and center are $e_j X_{L_i}^3$ and $(X_{L_i}^3 - 1, 0)$ respectively; the eccentricity of $L_i$ is provided by $e = \sqrt{\xi^2 + \eta^2}$. We note that the angular momentum $h$ of $L_i$ is constant. To represent the longitude of pericenter of the Lagrangian solutions as a function of $f_j$ we introduce the variables:

$$\tilde{\xi} := e \cos \varpi, \quad \tilde{\eta} := e \sin \varpi \quad (20)$$

which satisfy:

$$\begin{aligned} \tilde{\xi} &= \xi \cos(f_j + \varpi_j) + \eta \sin(f_j + \varpi_j) = e_j X_{L_i}^3 \cos \varpi_j + (X_{L_i}^3 - 1) \cos(f_j + \varpi_j) \\ \tilde{\eta} &= \xi \sin(f_j + \varpi_j) - \eta \cos(f_j + \varpi_j) = e_j X_{L_i}^3 \sin \varpi_j + (X_{L_i}^3 - 1) \sin(f_j + \varpi_j) \end{aligned} \quad (21)$$

i.e. $(\tilde{\xi}, \tilde{\eta})$ describes a circumference parameterized by $f_j$ whose radius and center are $|X_{L_i}^3 - 1|$ and $(e_j X_{L_i}^3 \cos \varpi_j, e_j X_{L_i}^3 \sin \varpi_j)$ respectively. We emphasize that $|X_{L_i}^3 - 1| \sim 3r_{Hill}$ where $r_{Hill} = \sqrt[3]{\frac{\varepsilon_j}{3}}$. The true anomaly $f$ and the longitude of perihelion $\varpi$ of $L_i$ are computed by inversion from $\xi, \eta$ and $\tilde{\xi}, \tilde{\eta}$. We note that<fnref n="3" />

$$\begin{aligned} &\text{if } 3r_{Hill} \gtrsim e_j \implies \varpi \in [0, 2\pi], f \text{ is bounded} \\ &\text{if } 3r_{Hill} \lesssim e_j \implies \varpi \text{ is bounded}, f \in [0, 2\pi]. \end{aligned} \quad (22)$$

Hence, for those planets fulfilling $3r_{Hill} \lesssim e_j$, the longitude of the perihelion of $L_1$ and $L_2$ is bounded on an interval whose length is

$$\Delta\varpi \sim 2 \arctan \frac{3r_{Hill}}{e_j}. \quad (23)$$

In Table 1 we write the mass, the semi-major axis, the eccentricity, the longitude of perihelion and $3r_{Hill}$ of Mercury, the Earth, and Mars. We note that for Mercury and Mars we have $3r_{Hill} \lesssim e_j$, hence the true anomaly of the Lagrangian points are such that $f \in [0, 2\pi]$ and $\varpi$ is bounded in some interval whose length is 0.11, and 0.30 respectively. In Figure 2 we plot the eccentricity, longitude of perihelion, and true anomaly associated to the Lagrangian solutions of the Sun–Mercury, Sun–Earth, and Sun–Mars ER3BPs computed through Eqs. 19 and 21. We remark that in the circular restricted three-body problem (CR3BP), the eccentricity and true anomaly of $L_1$ and $L_2$ are constants, while the longitude of perihelion varies linearly in $[0, 2\pi]$. Hence, one of the biggest difference between the CR3BP and the ER3BP is that the values of the longitude of perihelion associated to $L_1$ and $L_2$ can be bounded or not.

<footnote n="3">³We recall that all the dimensional quantities, such as the Hill radius and the semi-major axis, are expressed in these formulas in the rescaled units of measure so that the planet is at distance 1 from the Sun, referred as standard units.</footnote>
<table><caption>Table 1: Values of mass, semi-major axis, eccentricity, and longitude of perihelion used for the numerical experiments presented in this paper. On the sixth column we wrote the value of $3r_{Hill}$ in each Sun–planet standard unit (s.u.), On the seventh and eighth columns we wrote (in each Sun–planet standard unit) the minimum and maximum value of $\epsilon_j(1 + e_j \cos f_j)/(1 - e_j^2)$ with respect to the true anomaly $f_j$.</caption><thead><tr><th>Celestial body</th><th>Mass (kg)</th><th>a<sub>j</sub> (AU)</th><th>e<sub>j</sub></th><th>&omega;<sub>j</sub> (deg)</th><th>3r<sub>Hill</sub> (s.u.)</th><th>&epsilon;<sub>j</sub> (s.u.)</th><th>&epsilon;<sub>j</sub> (s.u.)</th></tr></thead><tbody><tr><td>Sun</td><td>1.9885 &times; 10<sup>30</sup></td><td>&mdash;</td><td>&mdash;</td><td>&mdash;</td><td>&mdash;</td><td>&mdash;</td><td>&mdash;</td></tr><tr><td>Mercury</td><td>3.302 &times; 10<sup>23</sup></td><td>0.387</td><td>2.06 &times; 10<sup>-1</sup></td><td>77</td><td>0.011</td><td>1.38 &times; 10<sup>-7</sup></td><td>2.09 &times; 10<sup>-7</sup></td></tr><tr><td>Earth</td><td>5.97237 &times; 10<sup>24</sup></td><td>1</td><td>1.67 &times; 10<sup>-2</sup></td><td>463</td><td>0.03</td><td>2.95 &times; 10<sup>-6</sup></td><td>3.05 &times; 10<sup>-6</sup></td></tr><tr><td>Mars</td><td>6.41712 &times; 10<sup>23</sup></td><td>1.52</td><td>9.336 &times; 10<sup>-2</sup></td><td>336</td><td>0.014</td><td>2.95 &times; 10<sup>-7</sup></td><td>3.56 &times; 10<sup>-7</sup></td></tr></tbody></table>

<figure><img src="image_3.png" /><figcaption>Figure 2: Evolution of the eccentricity, longitude of perihelion, and true anomaly of: $L_2$ in the Sun–Mercury ER3BP (top panels); $L_1$ (dashed line) and $L_2$ (continuous line) in the Sun–Earth ER3BP (middle panels); $L_1$ in the Sun–Mars ER3BP (bottom panels).</figcaption></figure>
**Computation of the orbital elements evolution for solutions in the manifold tubes.** To highlight the properties of the time evolution of the orbital elements for the solution in the manifold tubes we numerically integrate the equations of motion of the spacecraft using the ER3BP, with initial data computed through the Floquet-Birkhoff normal forms for different values of $f_j$ for the Sun-Mercury and Sun-Earth systems, and we calculate numerically the orbital elements using the swift package. For both systems we consider three orbits of the stable manifold tube of $L_2$ whose initial data are computed for $I_1 = 5 \times 10^{-7}$, $P_2 = 10^{-5}$, and $f_j = 1, 3.5$ and $6$; the orbits in the rotating-pulsating reference frame, their projection on the Oxy reference frame, and the orbital elements evolution are plotted in Figures 3 and 4. Accordingly to the previous discussion, we note that:

- when the spacecraft is close to the Lagrangian solution (referring to the labels of the Figures, when $\Delta t$ is close to 0) the orbital elements oscillate around those of the equilibrium;

- some time after the moment in which the spacecraft leaves the planet Hill sphere the orbital elements stabilize and the orbit is well approximated by a Keplerian one (up to an eventual close encounter with the planet). In the figures, we plotted in gray and pink the part of the orbit before such stabilization, in black and red the part of the orbit after such stabilization;

- the eccentricity and the longitude of perihelion stabilize quickly in all the cases. This implies that for those system characterized by $3r_{\text{Hill}} \lesssim e_j$, such as that of Mercury and Mars, the distribution of the longitude of the perihelion of the orbits in the manifold tubes, when they are far from the planet Hill sphere, is bounded.

- in the case of Mercury, the final values of $\varpi$ belong to an interval of amplitude comparable to the value $\Delta\varpi \sim 0.11$ provided by Eq. (23).

**Relevance for Hohmann transfers.** From a space flight dynamics point of view, this implies that if the orbits in the families of the manifold tubes of two planets $P_h, P_k$ have both bounded distribution of the longitudes of the perihelion, and if the intervals of these distributions do not match, it is not possible to perform an Hohmann transfer between the unstable manifold tube originating at a collinear point of planet $P_h$ to the stable manifold tube originating at a collinear point of the other planet $P_k$. In fact, the Hohmann transfer requires that the osculating ellipses of the departing and arrival points are co-axial. In such a case, other types of transfers may be considered, in the context of the ER3BP. If instead the two intervals partially overlap, Hohmann transfers may be found with suitable synchronization of the maneuvers in the context of the ER3BP. We notice that since for the Earth we have $3r_{\text{Hill}} \gtrsim e_h$, then the longitudes of the perihelion of the manifold tubes originating at $L_1$ and $L_2$ are distributed in $[0, 2\pi]$, and it is possible to determine an Hohmann transfer pushing the spacecraft from the unstable manifold tubes of the Earth to the stable manifold tubes of Mercury or Mars. In Section 3 we provide some sample computations.

**A test of the threshold effect.** We test the validity of the picture provided above for a fictitious planet having eccentricity $e_j$ close to the threshold value $\bar{e}_j \sim 3r_{\text{Hill}}$. We use as a model problem the gravitational parameter of the Earth, and we consider two values of $e_j$ slightly smaller or larger than $\bar{e}_j$, and we numerically compute backward in time the solutions of the ER3BP with initial data $X = X_{L_2} + 10^{-5}$, $Y = 0$, $P_X = 0$, $P_Y = X_{L_2}$ for different initial true anomalies of the planet, when eccentricities are $e_{\text{min}} = \bar{e}_j - \Delta e$, and $e_{\text{min}} = \bar{e}_j + \Delta e$, with $\Delta e = 5 \times 10^{-3}$; the orbital elements evolution of the numerically computed solutions are shown in Fig. 5. From the figure we note that for values of the planet eccentricity slightly smaller or larger than the threshold we obtain a complete different distribution of the longitude of perihelion associated to the manifold tube of the collinear Lagrangian points; in the first case ($e_{\text{min}} \sim 2.5 \times 10^{-2}$) the distribution is spread in $[0, 2\pi]$,
<figure><img src="image_4.png" /><figcaption>Figure 3: Representation of three orbits of the Sun–Mercury $L_2$ stable manifold tube in the rotating-pulsating and *Oxy* reference frame (top left and top right panel respectively), and the evolution of the angular momentum, eccentricity, and longitude of perihelion (bottom panels). The orbits were computed by integrating the initial data obtained through the Floquet-Birkhoff normal forms with $I_1 = 5 \times 10^{-7}$, $\phi_1 = 0$, $Q_2 = 0$, $P_2 = 10^{-5}$ and for different values of the initial true anomaly of Mercury; we choose $f_j = 1, 3.5, 6$. The position of the Sun and the initial data in the *Oxy* reference frame are plotted in the second panel through orange and black filled circles respectively. We plot in gray and pink the part of the orbits that are not in the Keplerian phase (see the previous discussion of the phases (i), (ii), (iii) for a solution in the manifolds tubes) and in black and red the part of the orbits that are in the Keplerian one.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Figure 4: Representation of three orbits of the Sun-Earth *L*₂ stable manifold tube in the rotating-pulsating and *Oxy* reference frame (top left and top right panel respectively), and the evolution of the angular momentum, eccentricity, and longitude of perihelion (bottom panels). The orbits were computed by integrating the initial data obtained through the Floquet-Birkhoff normal forms with *I*₁ = 5 × 10⁻⁷, *φ*₁ = 0, *Q*₂ = 0, *P*₂ = 10⁻⁵ and for different values of the initial true anomaly of Mercury; we choose *f*ⱼ = 1, 3.5, 6. The position of the Sun and the initial data in the *Oxy* reference frame are plotted in the second panel through orange and black filled circles respectively. We plot in gray and pink the part of the orbits that are not Keplerian, and in black and red the part of the orbits that are Keplerian. The dashed ellipse in the top right panel represents *L*₂ position in the *Oxy* reference frame.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Figure 5: Orbital elements evolution for different orbits in two ER3BPs where the planet $P_j$ has the same gravitational parameter of the Earth and eccentricities $e_{\min} = \bar{e}_j - \Delta e$ (continuous line) and $e_{\max} = \bar{e}_j + \Delta e$ (dashed line) with $\Delta e = 5 \times 10^{-3}$.</figcaption></figure>

in the second case ($e_{\max} \sim 3.5 \times 10^{-2}$) the distribution is bounded on an interval whose length is $\sim 1.16$. This interval is comparable with the interval $\Delta\omega \sim 1.42$ provided by Eq. (23).

**Time evolution of the angular momentum.** From the Figures 3 and 4 we notice that the final values for the angular momentum $h$ for all the orbits are very close (as a matter of fact, much closer than the final values of the eccentricity). In essence, this is a consequence of the fact that already on the Lagrangian solutions the eccentricity has a certain oscillation, while the angular momentum is constant. We therefore compute the time evolution of $h$ as the solution leaves a neighborhood of the Lagrangian one.

By representing the angular momentum and its time derivative using the rotating-pulsating canonical variables, we have:

$$h = xv_y - yv_x = \sqrt{1 - e_j^2} (P_Y X - P_X Y),$$

as well as:

$$\frac{dh}{dt} = \frac{dh}{df_j} \frac{df_j}{dt} = \{h, \mathcal{H}\} \frac{df_j}{dt} = \varepsilon_j \frac{1 + e_j \cos f_j}{1 - e_j^2} Y \left( 1 - \frac{1}{((X-1)^2 + Y^2)^{3/2}} \right). \quad (24)$$

The dependence of $\dot{h}$ on the orbital parameters of the planet appears as a factor multiplying the function:

$$w(X, Y) = Y \left( 1 - \frac{1}{((X-1)^2 + Y^2)^{3/2}} \right)$$

which depends only on the coordinates $X, Y$ of the rotating-pulsating orbital plane. With evidence, since $\dot{h}$ is proportional to $\varepsilon_j$, far from the planet $|\dot{h}|$ is small; but since the function $w(X, Y)$ is singular for $(X, Y) = (1, 0)$, close to the planet $\dot{h}$ is divergent; we also notice that we have $\dot{h} = 0$ for $Y = 0$ (so that we have $\dot{h} = 0$ on the Lagrangian solutions $L_1, L_2$), and for $(X - 1)^2 + Y^2 = 1$. We therefore check the magnitude of the values of $|\dot{h}|$ for $(X, Y)$ in the stable of unstable manifold tubes of interest.

In Fig. 6 we represent the values of $w(X, Y)$ using a color scale, and we plot the projection of some unstable manifold tubes of the Sun-Earth system (computed for $I_1 = 5 \times 10^{-7}$ and $f_3 = 0$) in the vicinity of $L_1$ and $L_2$; in the case shown by the figure, the maximum value of $|w(X, Y)|$ in the region visited by the manifold tubes is about $1.3 \times 10^3$; since the term $(1 + e_j \cos f_j)\varepsilon_j/(1 - e_j^2)$ is of order $3 \times 10^{-6}$ in the Sun-Earth case (the reader is referred to the seventh and eighth columns of Table 1), we obtain $|\dot{h}| \sim 4 \times 10^{-3}$. Moreover, we note that the manifold tubes are going away from the
regions characterized by the highest values of $|w(X, Y)|$ (these regions are colored in yellow). For such a reason the values of the angular momentum of the Lagrangian solutions and of the solutions $(X(f_j), Y(f_j), P_X(f_j), P_Y(f_j))$ in the branches of the manifold tubes which go far away from the planet remain very close up to an eventual additional close encounter with the planet. In Figure 7 we

<figure><img src="image_7.png" /><figcaption>Figure 6: Representation of the values of $|w(X, Y)|$ using a logarithmic color scale restricted to the interval $[10^{-1}, 2 \times 10^4]$. We plot also the projection of the unstable manifold tubes (white curves) of the Sun-Earth system whose initial data are computed via the Floquet-Birkhoff normal forms at both the Lagrangian points, for $I_1 = 5 \times 10^{-7}$, $f_j = 0$, $(Q_2, P_2) = (10^{-6}, 0)$ (for $L_1$), or $(Q_2, P_2) = (-10^{-6}, 0)$ (for $L_2$), and $\phi_1 \in [0, 2\pi]$.</figcaption></figure>

compare the range of the values of the angular momentum with the range of the orbital parameters $a, e$ in for a large set of orbits in the families $K_{I_1}^u, K_{I_1}^s$ of the Lagrangian points $L_1, L_2$ of planets Earth, Mars and Mercury which are of interest for interplanetary transfers from the Earth to Mars or Mercury (i.e. we consider $K_{I_1}^s$ of the Sun-Mercury $L_2$ and of the Sun-Mars $L_1$, and $K_{I_1}^u$ for both Sun-Earth $L_1, L_2$). The represented values have been obtained from numerical integrations of the ER3BP with initial conditions in the manifold tubes, computed using the Floquet-Birkhoff normal forms of order $N=6$, by fixing the values of $I_1$ and $Q_2$ (if we are considering the unstable manifold), or $P_2$ (if we are considering the stable manifold). On this occasion we generated the initial data by changing $\phi_1 \in [0, 2\pi]$ as well as the initial true anomaly of the planet $f_j \in [0, 2\pi]$; in such a way we are considering a large sample of the points of the planar torus $\mathcal{T}_{I_1}$. We generate initial data for two values of $I_1 = 5 \times 10^{-9}$ and $I_1 = 5 \times 10^{-7}$, using a uniform sampling of $\phi_1 \in [0, 2\pi]$ and $f_j \in [0, 2\pi]$. The solutions on the unstable manifold tube are numerically integrated for positive times, and the integration is stopped when the spacecraft is at the aphelion and is far from the planet; the solutions on the stable manifold tube are numerically integrated for negative times, and the integration is stopped when the spacecraft is at the perihelion and is far from the planet. The panels of 7 represent the values of the eccentricity and longitude of perihelion (vs the value of the angular momentum) of the spacecraft at these stopping points. We note that for the smaller value of $I_1$ we obtain a smaller interval for the angular momentum values (because greater amplitudes of the invariant torus represent larger oscillations around the equilibrium solutions), while the intervals for the eccentricity and the longitude of perihelion do not change so much by changing $I_1$.

## 3 Application to interplanetary missions

The previous theory of computation of the families $K_{I_1}^u, K_{I_1}^s$ of Keplerian orbits in the manifold tubes of the Lagrangian solutions $L_1, L_2$ of the ER3BP has been specifically formulated to design
<figure><img src="image_8.png" /><figcaption>Figure 7: Values of the angular momentum, eccentricity, and argument of perihelion at the stopping points for a large set of orbits in the family associated to the stable manifold tubes of the Sun–Mercury $L_2$ point (top panels), unstable manifold tubes of the Sun–Earth $L_1$ and $L_2$ points (middle panels), and the stable manifold tubes of the Sun–Mars $L_1$ point (below panels). The violet and green points refer the orbital elements of the tube manifolds computed with $I_1 = 5 \times 10^{-7}$ and $I_1 = 5 \times 10^{-9}$ respectively.</figcaption></figure>
interplanetary transfers. We illustrate an application using the simplest transfer maneuvers, i.e. the Hohmann transfers, between the orbits in the sets $K_{I_1}^u, K_{I_1}^s$ of the Lagrangian solutions $L_1$ or $L_2$ of the ER3BPs defined by two planets $P_h, P_k$: from the unstable manifolds of $L_1, L_2$ of the Sun–Earth ER3BP we transfer to the stable manifold of $L_2, L_1$ of the Sun–Earth and Sun–Mars ER3BPs respectively (we do not consider transfers from the Earth to Venus only because both planets have a small value of the orbital eccentricity, and the transfers can be conveniently designed using the two related CR3BPs). The initial conditions and the maneuvers of the transfers are computed using the ER3BPs, but we also check the effectiveness of the method with a numerical integration of the BER4BP, accounting of the gravitational interactions of the spacecraft with both planets $P_h, P_k$ during the full time interval needed for the spacecraft to move from a neighborhood of the Lagrangian solution of $P_h$ to a neighborhood of the Lagrangian solution of $P_k$.

The Hohmann transfer consists in two impulsive orbital maneuvers at two points A and B that are the aphelion and perihelion (or perihelion and aphelion) respectively of two osculating coplanar and co-axial ellipses in the families $K_{I_1}^u, K_{I_1}^s$.

We denote by A the aphelion or perihelion of the Keplerian orbit where the first maneuver is performed, and by B the perihelion or aphelion of the Keplerian orbit where the spacecraft is injected after the second maneuver. The possible values for the angular momentum, eccentricity and longitude of perihelion of A and B that can be used for the different transfers have been plotted in Fig. 7; we remark that, as explained in the previous Section, the possible longitudes of perihelion of the points A, B of the Sun–Mercury and Sun–Mars systems are distributed in a bounded interval $[\varpi_{\min}, \varpi_{\max}]$, while they spread in $[0, 2\pi]$ in the Sun–Earth case.

**Transfers in the two decoupled ER3BPs.** We denote by $S_A$ and $S_B$ the sets containing the angular momentum, eccentricity, longitude of perihelion, and the value of the true anomaly of the planet associated to the stopping points in the Sun–Earth unstable manifold, and in the Sun–Mercury (or Sun–Mars) stable manifold respectively. Each element in $S_A$ represents the aphelion of the Keplerian orbit associated to the unstable manifold tube in the Sun–Earth system computed by integrating in the future the ER3BP with initial data $r_{in}$ in the Oxy reference frame, and identifies the time-span $\Delta t_A$ needed by the spacecraft to depart from $r_{in}$ and reach A. Each element in $S_B$ represents the perihelion of the Keplerian orbit associated to the stable manifold tube in the Sun–Mercury or Sun–Mars system, computed by integrating in the past the ER3BP with initial data $r_{fin}$, and identifies the time-span $\Delta t_B$ needed by the spacecraft to depart from B and reach $r_{fin}$. Since the longitudes of perihelion associated both the manifold tube of the Sun–Mercury and Sun–Mars systems are distributed in a bounded interval $[\varpi_{\min}, \varpi_{\max}]$, while they are spread in $[0, 2\pi]$ in the Sun–Earth case (the reader is referred to Sec. 2 for further details), we decided to compute the optimal Hohmann transfer between the points defined by $S_A$ and $S_B$ by: i) selecting only the points $A' \in S_A$ characterized by longitudes of perihelion in the interval $[\varpi_{\min}, \varpi_{\max}]$; ii) for each point $A'$, selecting the point $B'' \in S_B$ whose longitude of perihelion is the closer one in the sampling of the set $S_A$: then we consider the point $B'$ having the same angular momentum and eccentricity of $B''$, the same longitude of perihelion of $A'$, and whose time-span to arrive at the vector state in the local stable manifold tube is $\Delta t_{B''}$; iii) once we determined all the possible pair $(A', B')$, we consider only the pair $(A, B)$ characterized by the minimum value of the Hohmann delta-v. In the first and third panel of Fig. 8 we plot the total time-span and the corresponding delta-v needed for the whole transfer between the Earth and Mercury (top panels) and the Earth and Mars (below panels) for all the pairs $(A', B')$ found in $S_A$ and $S_B$; for both transfers, we consider the sets $S_A$ and $S_B$ defined by $I_1 = 5 \times 10^{-7}$. In the second and fourth panel of Fig. 8 we plot the initial true anomalies at epoch $t=0$ with the corresponding delta-v needed for the transfer of the planets found for the Earth–Mercury and the Earth–Mars transfer respectively. We note from Fig. 8 that the difference between the maximum and the minimum delta-v is small (about 0.32
<figure><img src="image_9.png" /><figcaption>Figure 8: Top panels: values of delta-v found for the transfer between $\mathcal{K}_{I_1}^u$ of the Sun–Earth $L_1$ and $\mathcal{K}_{I_1}^s$ of the Sun–Mercury $L_2$ with $I_1 = 5 \times 10^{-7}$. Below panels: Values of delta-v found for the transfer between $\mathcal{K}_{I_1}^u$ of the Sun–Earth $L_2$ and $\mathcal{K}_{I_1}^s$ of the Sun–Mars $L_1$ with $I_1 = 5 \times 10^{-7}$. On the left of both the top and below panels, the total time interval needed for the whole transfer and the corresponding value of delta-v are shown; on the right, the values of delta-v for different initial true anomalies at epoch $t = 0$ of the planets at the beginning of the whole transfer. The cyan points indicate the values of the initial true anomalies corresponding to the minimum value of delta-v (about 15.543 km/s for the Earth–Mercury case, and 4.488 km/s for the Earth–Mars case).</figcaption></figure>
km/s and 0.12 km/s for the Earth–Mercury and the Earth–Mars transfer). In Table 2 the minimum values of delta-v in the Earth–Mercury transfers obtained by changing the value of $I_1$ are written. Once we determine the pair (A, B), we obtain the initial data $r_{in}$, i.e. the initial position vector

<table><caption>Table 2: Values of the minimum delta-v (in km/s) needed by an Hohmann transfer between $K_{I_1,E}^u$ of the Sun–Earth $L_1$ and $K_{I_1,M_e}^s$ for different values of $I_{1,E}$ and $I_{1,M_e}$, whose values are written in the first column and first row respectively.</caption><thead><tr><th>Sun–Mercury<br>Sun–Earth</th><th>5 × 10<sup>-7</sup></th><th>5 × 10<sup>-9</sup></th></tr></thead><tbody><tr><th>5 × 10<sup>-7</sup></th><td>15.543</td><td>15.581</td></tr><tr><th>5 × 10<sup>-9</sup></th><td>15.566</td><td>15.603</td></tr></tbody></table>

in the *Oxy* reference frame of the transfer belonging to the unstable Sun–Earth manifold tube, the time-span $\Delta t_A$ needed by the spacecraft to reach $r_A$ (the vector position in *Oxy* identified by A), the time-span $\Delta t_{tr}$ needed to transfer the spacecraft from $r_A$ to $r_B$ (the vector position in *Oxy* identified by B) through the Hohmann transfer, and the time-span $\Delta t_B$ needed to reach the final vector state $r_{fin}$. For the integration of the transfer we consider the equations of motion (2) where: i) from the starting point $r_{in}$ to $r_A$ we consider only the gravitational effects due to the Sun and the Earth; ii) from $r_A$ to $r_B$ we consider only the gravitational effects due to the Sun; iii) from $r_B$ to $r_{fin}$ we consider only the gravitational effects due to the Sun and Mercury (or Mars).

In Fig. 9 we plot the transfer between the Earth and Mercury; the manifold tubes in both Sun–Earth and Sun–Mercury systems were computed by defining $I_1 = 5 \times 10^{-7}$. In the top left panel, the transfer is shown in the *Oxy* reference frame, where the brown filled circle represent the location of Mercury at the end of the transfer, the black and cyan points represent $r_{in}$, $r_A$ and $r_B$. The panels below show the projection of the whole transfer in both the Sun–Earth and Sun–Mercury rotating-pulsating reference frame.

**Transfers in the BER4BP.** Since the vector positions $r_{in}$, $r_A$, $r_B$ and $r_{fin}$ are sufficiently far from Mercury (or Mars) and the Earth respectively, and also the gravitational parameters of the planets we are considering are small, when we integrate the equations of motion by considering the gravitational effects of the Sun, the Earth, and Mercury (or Mars) we expect that only small changes appear. For the computation of the orbital transfer, we decided to: i) integrate the Eqs. (2) considering both the Earth and Mercury (or Mars) with initial datum $r_{in}$, and stop the integration when the time-span is $\Delta t_A$. Let $r_A$ be the vector position of the stopping point; ii) integrate Eqs. (2) by using as initial vector position $r_A$ and as initial velocity vector the same Cartesian velocity vector found through Hohmann method for the transfer from A to B in the decoupled ER3BPs. The time-span of integration is $\Delta t_{tr}$. Let $r_B$ be the vector position of the stopping point; iii) integrate the Eqs. (2) with initial vector position $r_B$ and as velocity vector the same Cartesian velocity vector of the point B found in the arrival Sun–planet ER3BP. The integration time interval is $\Delta t_B$.

In the first panel of Fig. 10 we show the transfer between the Earth and Mercury using the BER4BP; both the manifold tubes of the Sun–Earth $L_1$ point and the Sun–Mercury $L_2$ points were computed using the Floquet-Birkhoff normal form with $I_1 = 5 \times 10^{-7}$. The transfer is almost equal to the transfer found using the ER3BPs (the reader is referred to the first panel of Fig. 9). This is because both Mercury and the Earth does not perturb so much the motion of the spacecraft when the latter is on the manifold of the Earth and Mercury respectively (the reader is referred to the right panel in Fig. 10 where the distance between the spacecraft and the Earth, and the spacecraft and Mercury during the whole transfer are plotted).
<figure><img src="image_10.png" /><figcaption>Figure 9: Transfer in the decoupled ER3BPs between the unstable manifold tube of the Sun–Earth $L_1$ (red orbit) and the stable manifold tube of the Sun–Mercury $L_2$ (blue orbit) in the *Oxy* reference frame (top left panel); for both system we compute the manifold tubes by imposing $I_1 = 5 \times 10^{-7}$ in the Floquet-Birkhoff normal forms. The brown filled circle represents the location of Mercury at the end of the transfer, while the black and cyan points denote the location of $r_{in}$, $r_A$, and $r_B$. The top right panel displays the distance between the spacecraft and the Earth (continuum line) and Mercury (dashed line) in the *Oxy* reference frame. The projection of the transfer in the Sun–Earth and Sun–Mercury rotating-pulsating reference frame is shown in the last two panels. For the integration we use the Sun–Earth standard unit of measure, i.e. the gravitational parameters of the Earth and the Sun are respectively $\varepsilon_3 = M_3/(M_0 + M_3)$, $\mu = 1 - \varepsilon_3$, the semi-major axis of the Earth is equal to one, while the gravitational parameters of Mercury and Mars are $\varepsilon_j = M_j/(M_0 + M_3)$ with $j = 1, 4$ ($M_0, M_1, M_3$ and $M_4$ are the masses of the Sun, Mercury, the Earth and Mars respectively; their values are represented in the second column of Table 1).</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Figure 10: Left panel: transfer between the Earth and Mercury in the BER4BP (for the integration we use the Sun-Earth standard unit). Right panel: distance between the spacecraft and the Earth (continuous line), and the spacecraft and Mercury (dashed line) in the *Oxy* reference frame during the whole transfer.</figcaption></figure>

## Acknowledgments

M. Guzzo and R.I. Paez acknowledge the project MIUR-PRIN 20178CJA2AB ‘New Frontiers of Celestial Mechanics: theory and applications’.

## References

E. Belbruno, M. Gidea, and F. Topputo. Weak stability boundary and invariant manifolds. *SIAM Journal on Applied Dynamical Systems*, 9:1061–1089, 2010.

M. Ceccaroni, A. Celletti, and G. Pucacco. Halo orbits around the collinear points of the restricted three-body problem. *Physica D*, 317:28–42, 2016.

H.D. Curtis. *Orbital mechanics for Engineering students*. Elsevier, 2014.

R.W. Farquhar, D.W. Dunham, Y. Guo, and J.V. McAdams. Utilization of libration points for human exploration in the sun-earth-moon system and beyond. *Acta Astronautica*, 55:687–700, 2004.

G. Gómez, À. Jorba, J. Masdemont, and C. Simó. *Dynamics and mission design near libration point orbits* Advanced Methods for Collinear Points, volume 3. Singapore: World Scientific, 2000.

G. Gómez, W.S. Koon, M.W. Lo, J.E. Marsden, J. Masdemont, and S.D. Ross. Connecting orbits and invariant manifolds in the spatial restricted three-body problem. *Nonlinearity*, 17:1571–1606, 2004.

M. Guzzo and E. Lega. On the identification of multiple close encounters in the planar circular restricted three-body problem. *Monthly Notices of the Royal Astronomical Society*, 428:2688, 2013.

M. Guzzo and E. Lega. A study of the past dynamics of comet 67p/churyumov gerasimenko with fast lyapunov indicators. *Astronomy and Astrophysics*, 579:A79, 2015.
M. Guzzo and E. Lega. Scenarios for the dynamics of comet 67p over the past 500 kyr. *Monthly Notices of the Royal Astronomical Society*, 469:S321, 2017.

M. Guzzo and E. Lega. Geometric chaos indicators and computations of the spherical hypertube manifolds of the spatial circular restricted three-body problem. *Physica D*, 373:38, 2018.

K.C. Howell and M. Kakoi. Transfers between the earth-moon and sun-earth systems using manifolds and transit orbits. *Acta Astronautica*, 59:367-380, 2006.

K.C. Howell, M. Beckman, C. Patterson, and D. Folta. Representations of invariant manifolds for applications in three-body systems. *The Journal of the Astronautical Sciences*, 54:69-93, 2006.

À. Jorba and J. Masdemont. Dynamics in the center manifold of the restricted three-body problem. *Physica D*, 132:189-213, 1999.

W.S. Koon, M.W. Lo, J.E. Marsden, and S.D. Ross. *Dynamical systems, the three-body problem and space mission design*. Marsden Books, 2008.

R.I. Paez and M. Guzzo. A study of temporary captures and collisions in the circular restricted three-body problem with normalizations of the levi-civita hamiltonian. *International Journal of Non-Linear Mechanics*, 120:103417, 2020.

R.I. Paez and M. Guzzo. Transits close to the lagrangian solutions $l_1, l_2$ in the elliptic restricted three-body problem. *Nonlinearity*, 34:6417, 2021.

J. Palacián and P. Yanguas. Invariant manifolds of spatial restricted three-body problems: the lunar case *New Advances in Celestial Mechanics and Hamiltonian Systems*, 2004.

J. Parker and R.L. Anderson. *Low-energy lunar trajectory design*. Pasadena, CA: Jes Propulsion Laboratory, California Institute of Technology, 2013.

P. Pergola, K. Geurts, C. Casaregola, and M. Andrenucci. Earth-mars halo to halo low-thrust manifold transfer. *Celestial Mechanics and Dynamical Astronomy*, 105:19-32, 2009.

G. Pucacco. Structure of the centre manifold of the $l_1, l_2$ collinear libration points in the restricted three-body problem. *Celestial Mechanics and Dynamical Astronomy*, 131, 2019.

Y. Ren, J.J. Masdemont, G. Gòmez, and E. Fantino. Two mechanisms of natural transport in the solar system. *Communication in Nonlinear Science and Numerical Simulation*, 17:844-853, 2012.

C. Simó. *Hamiltonian Systems with Three or more degrees of freedom*, volume 533. NATO Advanced Science Institutes Series C. Mathematical and Physics Sciences, 1999.

F. Topputo, M. Vasile, and F. Bernelli-Zazzera. Low energy interplanetary transfers exploiting invariant manifolds of the restricted three-body problem. *The Journal of the Astronautical Sciences*, 53:353-372, 2005.

A. Zanzottera, G. Mingotti, R. Castelli, and M. Dellnitz. Intersecting invariant manifolds in spatial restricted three-body problems: design and optimization of earth-to-halo transfers in the sun-earth-moon scenario. *Communications in Nonlinear Science and Numerical Simulations*, 17:832, 2012.