Characterization of Low-Energy Quasiperiodic Orbits in the Elliptic Restricted 4-Body Problem with Orbital Resonance

In this work, we investigate the behavior of low-energy trajectories in the dynamical framework of the spatial elliptic restricted 4-body problem, developed using the Hamiltonian formalism. Introducing canonical transformations, the Hamiltonian function in the neighborhood of the collinear libration point L1 (or L2), can be expressed as a sum of three second order local integrals of motion, which provide a compact topological description of low-energy transits, captures and quasiperiodic libration point orbits, plus higher order terms that represent perturbations. The problem of small denominators is then applied to the order three of the transformed Hamiltonian function, to identify the effects of orbital resonance of the primaries onto quasiperiodic orbits. Stationary solutions for these resonant terms are determined, corresponding to quasiperiodic orbits existing in the presence of orbital resonance. The proposed model is applied to the Jupiter-Europa-Io system, determining quasiperiodic orbits in the surrounding of Jupiter-Europa L1 considering the 2:1 orbital resonance between Europa and Io.


Introduction
Low-energy trajectories have long been used in space exploration, starting from missions to the libration points of the Sun-Earth system [1][2][3][4][5] and including lunar transfers passing close to the intermediate Earth-Moon libration point L 1 [6][7][8], named internal transfers, or extending beyond the orbit of the Moon to take advantage of the Sun gravitational attraction [9,10], named external transfers.
The spatial Circular Restricted 3-Body Problem (CR3BP) [11] represents an effective model for the preliminary analysis of both missions to the libration points [12][13][14][15] and low-energy transfers [16][17][18]. As proved by Conley [19], in this dynamical framework the ultimate behavior of low-energy trajectories can be determined based on their phase space description in the neighborhood of the collinear libration points [20], a property which has been extensively used in mission oriented works, since it allows a rapid design of low-energy missions [21][22][23][24][25].
More accurate solutions can be generated based on more sophisticated models, such as the Elliptic Restricted 3-Body Problem (ER3BP) [26], in which the primaries move along elliptic orbits, or the Circular Restricted 4-Body Problem (CR4BP), in which the spacecraft dynamics evolves under the gravitational attraction of four primaries (i.e., the Sun-Earth-Moon system) [27]. Unfortunately, because no integrals of motion, such as the Jacobi constant for the CR3BP, are known for neither the ER3BP nor the CR4BP, a general characterization of low-energy trajectories in these dynamical frameworks has not been derived yet.
In fact, the phase space description derived for the CR3BP in the surrounding of the libration point L 1 (or L 2 ) can actually provide pivotal results for the design of missions also in the presence of non negligible, but limited, eccentricity in the motion of the primaries [28] and gravitational attraction from a fourth body [29]. This result has its foundation on a theorem by Conley and Easton [30], stating that the basic topological properties of the phase space flow of the CR3BP are persistent in the presence of perturbations, whose validity was verified also by means of numerical analyses [31,32].
In this paper, we present a novel model that allows characterizing low-energy trajectories for a system consisting of a central body and two natural satellites moving along elliptic orbits, therefore in the dynamical framework of the spatial Elliptic Restricted 4-Body Problem (ER4BP). Typically, this task is performed based on the results of the CR3BP, corrected and compensated to take into account the perturbations introduced by the eccentric motion of the primaries or by the gravitational attraction of a third primary. The model proposed in this paper provides results which do not require any further correction and, furthermore, allows determining new solutions corresponding to quasiperiodic orbits existing when the primaries are in orbital resonance.
The ER4BP is modeled using the Hamiltonian formalism, where canonical transformations of the state variables allow reducing the Hamiltonian function for the ER4BP to a form equivalent to that of the CR3BP in the neighborhood of L 1 plus higher order terms, thus the complexity of the problem is transferred in the definition of a coordinate transformation. Following this approach, the low-energy trajectories can be characterized in terms of the new phase space variables, using the same techniques developed for the CR3BP [33], then the corresponding physical (position and velocity) coordinates can be derived applying the inverse coordinate transformation.
The effect of orbital resonance between two of the primaries is finally examined. Under this circumstance, any attempt of introducing a canonical transformation to absorb the higher order terms defaults, because of the small denominators problem [34], and the Hamiltonian function for the ER4BP, is transformed to that of the CR3BP plus an additional resonant term.
Specific phase space states in the transformed coordinates correspond to equilibrium conditions for the resonant term and these are actually periodic solutions in physical coordinates. The problem is here applied to the specific case of the Jupiter-Europa-Io system, in which the two moons are in a 2:1 resonance, to determine quasiperiodic trajectories in the neighborhood of Jupiter-Europa L 1 . Here, a satellite can be deployed to perform, or support, space science missions to the Jovian moon, which has attracted the interest of space agencies because of its geodynamically active icy crust, which makes Europa a major candidate for the search of life in the solar system [35,36].
The paper is organized as it follows, in Section 2 a model for the ER4BP is developed in space and velocity coordinates and converted to Hamiltonian variables. In Section 3 the Hamiltonian function for the ER4BP is set in a form equivalent to that of the CR3BP by means of canonical transformations. In Section 4 the resonant terms, depending on the true anomalies of the primaries, are identified, by means of small denominators problem, and the corresponding quasiperiodic solutions are determined. An application to the Jupiter-Europa-Io system is provided in Section 5, comparing the analytical solution obtained in Section 4 to the numerical results obtained from integrating the full nonlinear equations of motion associated to the ER4BP.

Dynamical Model
The dynamics of a restricted 4-body system, with masses m 1 , m 2 , m 3 >> m, is hereafter described under the following hypotheses: • The center of mass of the system is indicated as O; • The relative motion of m 2 with respect to m 1 describes a Keplerian orbit onto the plane Π with semimajor axis a and eccentricity e; • The relative motion of m 3 with respect to m 1 describes a Keplerian orbit onto the plane Π p with semimajor axis a p and eccentricity e p ; • The orbital plane Π p is tilted of an angle ε p with respect to Π; where the subscript p is used to explicit the parameters referring to the perturbing fourth body (m 3 ). This model is suitable to represent a multi-moon environment, which is typical for the outer planets, or gas giants, in the solar system. An inertial reference frame F i = Ξ ,Ĥ,Ẑ is introduced. F i is centered in O, withẐ orthogonal to Π,Ξ parallel to the line that connects m 1 and m 2 at the initial time t 0 , andĤ completing the rectangular reference frame. For the sake of simplicity, it is assumed that at t 0 the primary m 3 is at the apocenter of its orbit around m 1 and lies onto the Ξ ,Ẑ plane. It shall be noted that this hypothesis does not cause any loss of generality.
The dynamic equations of motion for m in the inertial reference frame are reported below where (Ξ i , H i , Z i ) indicate the coordinates of the i-th primary, and R i indicates the distance between m and the i-th primary. System (1) can be conveniently converted into rotating-pulsating coordinates, after introducing the reference frame F r = ξ ,η,ζ , which rotates rigidly with the primaries m 1 and m 2 . In particular, F r is parallel to F i at t 0 , with axisξ pointing from m 1 to m 2 , axisζ orthogonal to Π and axisη completing the orthogonal reference frame. As sketched in Figure 1 (the trajectory of the primaries in Figures 1 and 2 does not correspond to an actual system and is produced for the sake of representation, selecting: a = 6.4e + 5 km, a p = 4.2e + 5 km, e = 9.4e − 3, e p = 4.1e − 1, ε p = 43 deg, m 1 = 1.9e + 25 kg, m 2 = 8.9e + 24 kg, m 3 = 4.8e + 24 kg), in the rotating frame the primaries m 1 and m 2 oscillate along the axiŝ ξ. The trajectory of the primary m 3 is more complex and results from the combination of the periodic motion of m 3 around O and the periodic rotation of F r . Converting the inertial coordinates into rotating ones produces the following form for the dynamic equations of motion [37]  and θ is the true anomaly of m 2 .  Introducing the normalization factors for the mass M = and time T = R 2 GM , allows transforming the rotating coordinates into rotating-pulsating ones (x, y, z), and, finally, rearranging System (2) in a dimensionless form, with the dependence on time being replaced by that on θ [38] M and the apostrophe indicates the derivative with respect to θ. As shown in Figure 2, the position of the primaries m 1 and m 2 in rotating-pulsating coordinates corresponds to, respectively, [−µ, 0, 0] and [1 − µ, 0, 0] with µ = m 2 M . As in the rotating frame, the primary m 3 evolves along a three-dimensional trajectory which in this case depends also on the periodic motion of the rotating-pulsating reference frame along axis ξ.
The second order equations in System (3) are typically expressed as a set of first order equations for the position (x, y, z) and velocity v x , v y , v z coordinates The reader will notice that setting e = e p = µ 3 = 0 the Systems (1)-(4) reduce to the corresponding ones for the CR3BP. In Section 3 it will be shown how to express the dynamics of the ER4BP to a form equivalent to that of the CR3BP also when e, e p , and µ 3 are small but not negligible. To achieve this goal, the Hamiltonian formalism is introduced.

Hamiltonian Formalism
System (4) can be represented using the Hamiltonian formalism after deriving the Hamiltonian function for the ER4BP (see [26,39]) which depends on the variables q and p, where the bold style is used to indicate the vectors, named the positions and conjugate momenta The dynamic equations of motion are obtained by the definitions [40] q = ∂H ∂p and p = − ∂H ∂q , leading to Introducing Equation (6) into Equation (7) results in System (4), proving the equivalence between the Hamiltonian and the position space representations. Considering the non-polynomial terms of H F e, e p , the dependence on e, e p , and  (2) e, e p , µ 3 /r 3 (9) with where o (2) collects the terms of order higher than one in e, e p , and µ 3 /r 3 . For the sake of clarity, the long expressions for the partial derivatives of ρ i = r 2 i are reported in Appendix A. Collecting the polynomial terms of Equation (5) and F , results in the Hamiltonian of the CR3BP [11], hereafter indicated as Then, Equation (5) can be expressed in the compact form

Persistence of the Topological Properties
In Section 2, the Hamiltonian function of the ER4BP was expressed as the sum of the Hamiltonian for the CR3BP (H c ) and the terms depending on the orbital eccentricity of the primaries (e, e p ) and on the gravitational attraction of the third primary (µ 3 /r 3 ).
It is worth recalling now that for the CR3BP, the ultimate behavior of low-energy trajectories can be determined based on the topological properties of the linear flow in the phase space surrounding the collinear libration points L 1 (or L 2 ) [20], called the equilibrium region. According to a theorem by Conley and Easton, the properties of the flow in the equilibrium region are persistent also in the presence of perturbations [30], therefore for small e, e p and µ 3 /r 3 , it shall be possible to transform the Hamiltonian function of the ER4BP to a form equivalent to that of the CR3BP in the equilibrium region, hereafter indicated as H (2) c . The expression of H (2) c is derived hereafter. First, the origin of the system is translated to the libration point (equivalent solutions can be derived for the second libration point by simply replacing L 1 with L 2 ), setting Then the non polynomial terms of Equation (14) are rearranged as follows and expanded in power series up to order two in q, obtaining where L 1 is the x coordinate of the libration point, µ = m 2 m 1 +m 2 and Introducing Equations (18) and (19) into Equation (15) and excluding terms inessential in the Hamilton equations of motion leads to with In the following sections, two canonical transformations are introduced. The first one leads the second order Hamiltonian function for the ER4BP to a form equivalent to Equation (20), the second one allows expressing the resulting second order Hamiltonian function as the sum of three local integrals of motion.

Normal Forms for the Elliptic Restricted 4-Body Problem
A canonical transformation from the old to a new set of Hamiltonian variables T 1 : (q, p) → (Q, P) is designed by means of the generating function S, which verifies the conditionĤ = H(q, p, θ) + ∂S ∂θ (21) and links the old and new variables as The authors proved that selecting the generating function reported below (see Appendix B) The canonical transformation can absorb the second order terms in Equation (15) that depend on the perturbations (e, e p , µ 3 /r 3 ), hereafter explicited (please note that terms depending on the perturbation e p do not appear in Equation (23) because they are of a order higher than two) Then, the second order HamiltonianĤ (2) in the new variables reduces toĤ (2) and the full Hamiltonian is given bŷ c (Q, P) +Ĥ (n) Q, P, e, e p , µ 3 /r 3 (26) whereĤ (n) collects terms of order 3 or higher.
It is worth to summarize the result here. Applying the canonical transformation defined in Equation (21) with the generating function described by Equation (23), the second order Hamiltonian functionĤ (2) for the ER4BP results in Equation (25), whose form is equivalent to that of the second order Hamiltonian H (2) for the CR3BP, given by Equation (20). Because the topological properties of the phase space flow in the equilibrium region depend on the local form of the Hamiltonian function [20], the equivalence in the form ofĤ (2) and H (2) indicates that the topological properties of the CR3BP are equivalent to those of the ER4BP, substantiating the theorem by Conley and Easton [30].

Topological Characterization of Low-Energy Trajectories
A second canonical transformation T 2 : (Q, P) → (x, y), originally proposed by Siegel and Moser [41] and extensively used for the CR3BP [29,33,42], allows expressingĤ (2) as the sum of three terms where ρ, λ 1 , and λ 2 are, respectively, the real and the two imaginary eigenvalues of the linear system associated toĤ (2) and h indicates the constant value ofH (2) , named the energy level.
As proved by Moser [43], the three terms in Equation (27) are local integrals of motion, each one depending on a different couple of variables (x i , y i ). The two integrals associated with the imaginary eigenvalues λ i represent harmonic oscillators and characterize, respectively, the in-plane and out-of-plane oscillations, corresponding to the (x, y) and z components in position space coordinates.
The level surface ρx 1 y 1 describes the linear phase space flow in the neighborhood of the saddle point (x 1 = 0, y 1 = 0), as sketched in Figure 3. Low-energy trajectories can be characterized based on their state projection in (x 1 , y 1 ) and are classified as follows [20]: • Lissajous quasiperiodic orbits, characterized by x 1 = y 1 = 0, which evolve inside the equilibrium region; • Transit trajectories, corresponding to the hyperbolic segments x 1 y 1 < 0, which cross the equilibrium region twice, once towards m 1 and once towards m 2 , in a finite interval of time; • Bouncing trajectories, corresponding to the hyperbolic segments x 1 y 1 > 0, which never cross the equilibrium region; • Long-term ballistic captures, characterized by either x 1 → 0 or y 1 → 0, which cross the equilibrium region twice, over an indefinitely long interval of time.
Inasmuch, some relevant results proved that the CR3BP can be extended to the ER4BP, indicating that correlations exist between the osculating orbital elements of ballistic captures in the proximity of the primaries and their state (x, y) in the equilibrium region [33,44].
As detailed in the next section, the model proposed can provide further information on the dynamical evolution of quasiperiodic orbits under the effect of perturbations included in higher order terms. In particular, the impact of orbital resonance between the primaries is examined.

Identification of Resonant Terms
The focus now is on investigating the effects of residual perturbations, included in the higher order terms, onto quasiperiodic orbits (x 1 = y 1 = 0). In particular, the case of orbital resonance between the primaries is examined, which allows settingθ p = nθ, where n is an integer number.
The Hamiltonian function including the higher order terms can be easily derived, by applying the canonical transformation T 2 by Siegel and Moser toĤ (n) , producingH (n) , and adding this term to Equation (27), resulting in A new canonical transformation T 3 : (x, y) → (w, z), aiming at absorbing the order three terms in Equation (28), is now introduced. Operating as in Section 3, the canonical transformation is defined by means of a generating function S, which verifies H =H(x, y, θ) + ∂S ∂θ (29) and y = ∂S ∂x w = ∂S ∂z (30) Using symbolic algebra, the authors proved that setting the following expression for the generating function (see Appendix B) S(x, z, θ, g i ) = x 1 z 1 + x 2 z 2 + x 3 z 3 + g 1 x 2 2 + g 2 x 2 3 + g 3 z 2 2 + g 4 z 2 3 + g 5 z 2 + g 6 x 2 + Allows absorbing the order three terms of Equation (28), leading tō whereH (m) collects the terms of order higher than three. As expected from the small denominator problem [34], in case of resonance, the denominators of some coefficients in g i are equal to zero and the canonical transformation defaults. In particular, this occurs in g 5 and g 6 for n = λ 1 + and in g 7 and g 8 for n = λ 2 + , where is an arbitrarily small constant.
The residual order three terms of the Hamiltonian function are reported below where j = 2, 3 and b

Determination of Stationary Points
The model developed here allows determining periodic solutions existing because of orbital resonance. Before proceeding with their determination, the following property of the canonical transformation by Siegel and Moser shall be recalled [40,41] Introducing Equation (34) into Equation (33) and applying Euler's formula, to express the exponential terms as trigonometric functions, the expression ofH res j is transformed tõ Here the dependence on the true anomaly θ can be absorbed by introducing the canonical transformation T 4 : (ρ j , φ j ) → (ρ j ,φ j ) defined by the following generating function Corresponding to the change of coordinates reported below (the variables not reported are not transformed) Which transforms Equation (35) as follows Stationary points exist for Equation (38), corresponding to The reader shall note that the sequence of transformations from T 1 to T 4 , hereafter indicated as T 1,4 and producing [ρ, σ] = T 1,4 [q, p], depend on the commensurable θ and θ p . Applying the inverse of T 1,4 to any stationary solution defined by Equation (39) results in a quasiperiodic trajectory in [q, p] s and related position and velocity components.

Numerical Analysis
In the previous sections, we derived a model which allows designing quasiperiodic libration point orbits existing in case of resonance between the primaries. This model is here verified by means of numerical analysis on the Jupiter-Europa-Io system, whose parameters are reported in Table 1 (from the values in Table 1 it is possible to compute µ 3 /r 3 << 10 −3 ). This system is characterized by the 2:1 orbital resonance between Europa and Io, corresponding to n = 2.02.
The effects of the orbital resonance is investigated on quasiperiodic orbits in the neighborhood of Jupiter-Europa L 1 , characterized by λ 1 = 2.22 and λ 2 = 2.15. Since n = 2.02, the condition n − λ i = is verified. In particular, the case of an out of plane periodic orbit is presented below. The initial conditions, selected from Equation (39) for j = 3, arẽ φ 3,s = 0 andĤ 3 res = 10 −15 and correspond to x 1 = L 1 , x 2 = 0, x 3 = −3.11 × 10 −10 , y 1 = 0, y 2 = 0 and y 3 = 4.556 × 10 −3 . From this initial state, the full nonlinear equations of motion in System (4) are integrated using ode113 solver in Matlab, for a total time equal to 200 times the orbital period of Europa around Jupiter, which corresponds to about 708 Earth days. The time behavior of the oscillations is shown in Figure 4, where the distance is normalized by DU = R = a.  The solution is compared to the one obtained ignoring the order three terms, corresponding to the state x 1 = L 1 , x 2 = 0, x 3 = −10 −10 , y 1 = 0, y 2 = 0 and y 3 = 0. The time behaviour of the oscillations is represented in Figure 5, showing that the orbital perturbation of Io produces an increase of the amplitude up to ±2.9 × 10 −4 DU in the time interval considered. Differently, the amplitude of oscillations for the quasiperiodic orbit obtained from the stationary point never exceeds the value 10 −5 DU. A representation of the quasiperiodic trajectory in the rotating frame is shown in Figure 6. It can be observed that the vertical component oscillates as in Figure 4, while the in-plane components have a negligible drift towards the negative values, resulting from higher order perturbations. It shall be emphasized that over the 200 orbits, the drift of the in-plane components is four orders of magnitude lower than the amplitude of the vertical oscillations. To provide a better insight on the shape of the quasiperiodic orbit corresponding to the stationary point, the trajectory is represented in Figure 7 in the inertial reference frame (Ξ, H, Z).

Conclusions
In this manuscript we presented a model to extend the topological characterization of low-energy trajectories from the dynamical framework of the CR3BP to that of the ER4BP. This goal is achieved, defining a sequence of canonical transformations which convert the Hamiltonian function of the ER4BP, evaluated in the neighborhood of the libration point L 1 (or L 2 ) to a form equivalent to that of the CR3BP.
The higher order terms of the Hamiltonian function in the normal form, which depend on the true anomalies of the primaries, can be absorbed by a further canonical transformation. This transformation defaults in case of orbital resonance between the primaries, as indicated by the problem of small denominators.
For the terms of the Hamiltonian function that cannot be absorbed by the canonical transformation, in virtue of the orbital resonance, stationary points can be computed, which correspond to quasiperiodic orbits in the position space coordinates. The above mentioned model was verified, by means of numerical analyses, to investigate the behavior of quasiperiodic orbits for the Jupiter-Europa-Io system, characterized by the 1:2 resonance of Io and Europa. The full nonlinear equations of motion for the ER4BP were integrated for an initial condition corresponding to the vertical harmonic oscillator at Jupiter-Europa L 1 .
The results clearly showed that neglecting the terms of order three results in a trajectory in which the amplitude of the oscillations increases indefinitely in time, driven by the orbital resonance between the primaries. Instead, using the initial conditions associated to stationary points determined for the order three resonant Hamiltonian, reported in Equation (27), the amplitude of oscillations does not diverge and the resulting orbits are quasiperiodic.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. Terms of the Power Series Expansion
The expressions for the partial derivatives of ρ i , that appear for the first time in Equations (10)- (13), are reported below, in order to allow the reader to replicate the results presented. The derivatives that are not reported are equal to zero.