Issue 
A&A
Volume 638, June 2020



Article Number  A137  
Number of page(s)  10  
Section  Celestial mechanics and astrometry  
DOI  https://doi.org/10.1051/00046361/202037696  
Published online  26 June 2020 
Shortperiod effects of the planetary perturbations on the Sun–Earth Lagrangian point L_{3}
Planetary perturbations of the Sun–Earth L_{3}
^{1}
Università degli Studi di Padova, Centro di Ateneo di Studi e Attività Spaziali “Giuseppe Colombo” – CISAS, Via Venezia 15, Padova 35131, Italy
email: erica.scantamburlo@phd.unipd.it
^{2}
Università degli Studi di Padova, Dipartimento di Matematica “Tullio LeviCivita”, Via Trieste, 6335121 Padova, Italy
email: guzzo@math.unipd.it
Received:
10
February
2020
Accepted:
9
April
2020
Context. The Lagrangian point L_{3} of the Sun–Earth system, and its Lyapunov orbits, have been proposed to perform stationkeeping, although L_{3} is only rigorously defined for the extremely simplified model represented by the reduced Sun–Earth–spacecraft system. As in L_{3} the planetary perturbations (mainly from Jupiter and Venus) are stronger than Earth’s attraction, it is necessary to understand whether or not the dynamics close to L_{3} persist under such a strong perturbation, allowing for a definition of dynamical substitutes for models that are more realistic than the circular restricted threebody problem.
Aims. In this paper we address the problem of the existence of motions that remain close to L_{3} for a timespan which is relevant for space missions in a model of the Solar System compatible with the precision of JPL digital ephemerides.
Methods. First, we computed analytically the main shortperiod effects of planetary perturbations in a simplified model of the Solar System with the orbits of all the planets coplanar and circular. We then applied the Fast Lyapunov Indicator method in order to find dynamical substitutes that exist for timespans of hundreds of years in the model of the Solar System that is used to produce the modern ephemerides.
Results. We find that the original system is conjugate by a canonical transformation to an averaged system that has an equilibrium close to L_{3}: even if Venus and Jupiter each move the position of this equilibrium by about 218 and 176 km, respectively, in opposite directions, in the model where both the planets are included, their effects almost perfectly compensate for one another, leaving a displacement of about 40 km only. This equilibrium is then mapped in the original system to a quasiperiodic dynamical substitute; the contributions of each planet to the amplitude of this quasiperiodic libration around L_{3} do not compensate for one another, and sum to about 10 000 km. The Fast Lyapunov Indicator method allowed us to find orbits of any amplitude bigger than this one (up to 0.03 AU) for timespans of hundreds of years in the model of the Solar System that is used to produce the modern ephemerides.
Conclusions. Using a combination of the Hamiltonian averaging method with a new implementation of the Fast Lyapunov Indicator method we find orbits useful for astrodynamics originating at the Sun–Earth Lagrangian point L_{3} for a realistic model of the Solar System. In particular, this usage of the chaos indicator provides an innovative application of dynamical systems theory to astrodynamics, where the shortperiod perturbations represent a relevant part of the model.
Key words: celestial mechanics
© ESO 2020
1. Introduction
The dynamics close to the Lagrangian points L_{1} and L_{2} of the Sun–Earth circular restricted threebody problem (CR3BP) have been extensively discussed in recent decades in connection with spacemission design (the reader is referred to Martin et al. 2010; van Damme et al. 2010; Koon et al. 2006; Parker & Anderson 2013; Gómez et al. 2004, and references therein for further information). ISEE3, SOHO, WMAP, Genesis, Herschel and Planck are just few examples of missions in which the spacecraft reached the vicinity of the L_{1} and L_{2} Lagrangian points of the Sun–Earth system during its journey. While L_{3} is also potentially useful for astrodynamics (Gómez et al. 2001; Gómez & Mondelo 2001; Barrabés & Ollè 2006; Terra et al. 2014; Pàez & Efthymiopoulos 2015; Jorba & Nicolás 2020), no significant attention has been paid to the Sun–Earth L_{3} point even if its Lyapunov orbits have been proposed as an ideal place for spacecraft because of the possibility to obtain data for the Sun, the inner planets, and the main belt from a new perspective (Hou et al. 2007; Tantardini et al. 2010). The aim of this study is to find orbits close to the Sun–Earth L_{3} point that are useful for applications of astrodynamics. In particular, this refers to orbits that remain close to L_{3} within useful distances at least for several decades in a model of the Solar System whose ephemerides are comparable with those provided by the JPL computation service. Because the planetary perturbations on the Sun–Earth Lagrangian point L_{3} are stronger than Earth’s attraction itself, the CR3BP is not a realistic approximation of the dynamics close to L_{3}. The main differences between the real Solar System and the Sun–Earth circular restricted problem are caused by the planetary perturbations (mainly from Venus and Jupiter) and in the elliptic orbit of the Earth. To quantify the effects of planetary perturbations, we first consider a multicircular restricted problem (MCR8BP), where the eight planets of our Solar System (from Mercury to Neptune) perform circular motions around the Sun. Using Hamiltonian perturbation theory we find that the Hamiltonian of the MCR8BP is conjugate by a canonical transformation to an averaged Hamiltonian system that still has an equilibrium point close to L_{3} (within terms of second order in the planetary masses and fourth order in eccentricity). This equilibrium, which we denote by ℒ_{3}, is aligned in the Sun–Earth direction. The individual contribution of each planet to the displacement of ℒ_{3} from L_{3} is reported in the first column of Table 1. It is remarkable that while Venus and Jupiter each move the position of this equilibrium by about 218 and 176 km, respectively, in opposite directions, in the model where both planets are included, their effects almost perfectly compensate for one another, leaving a displacement of about 40 km only (the contributions of the other planets being smaller by about two orders of magnitude). The equilibrium ℒ_{3} is then mapped in the original system to a quasiperiodic dynamical substitute; the contributions of each planet to the amplitude of this quasiperiodic libration around L_{3}, which are reported in the second column of Table 1, are quite large and sum to about d = 10^{−4} AU ∼15 000 km. Since the linear stability of ℒ_{3} is the same as that of L_{3}, the averaged Hamiltonian, approximated within terms of second order in the planetary masses, has close to ℒ_{3} a dynamics similar to that of L_{3} in the CR3BP; this means that in the vicinity of ℒ_{3} there are periodic and quasiperiodic orbits. The image of these orbits in the original MCR8BP provides orbits librating around L_{3} of amplitudes larger than the amplitude of the dynamical substitute, which we denote hereafter by d. Since the MCR8BP is still a rough approximation of the real Solar System, we use fast Lyapunov indicators (FLI) to look for orbits of this kind in a model of the Solar System that is compatible with the precision of JPL digital ephemerides. In the last decade, the FLI have been used to find orbits related to the dynamics close to the Lagrangian points of the restricted threebody problem (Lega et al. 2011; Guzzo & Lega 2013, 2014, 2018) and to close encounters in cometary dynamics (Guzzo & Lega 2013, 2017). The method provides orbits that remain, within a given timespan, in a neighborhood ℬ of the dynamical substitute obtained from an approximated model. In order to construct the neighborhood ℬ, we must also consider the differences between the real Solar System and the Sun–Earth CR3BP caused by the elliptic orbit of the Earth. The CR3BP is traditionally compared to the elliptic restricted threebody problem (ER3BP) by introducing a rotatingpulsating Cartesian reference frame in which the primary and the secondary bodies remain in fixed locations of the axis (see e.g., Szebehely 1967). More precisely however, in the rotatingpulsating reference frame the ER3BP still has the Lagrange equilibria in the same positions as the circular problem. We therefore proceed as follows: for any given libration amplitude α > d, we first compute a Lyapunov orbit of the CR3BP of amplitude α, hereafter referred to as the planar Lyapunov orbit. We use this orbit as an input of the FLI method to find orbits of the ER3BP that, in the pulsating frame, remain close to the previous Lyapunov orbit for up to 200 years. We then consider the motions of the Sun and the Earth in the realistic model of the Solar System, and we introduce a new rotatingpulsating reference frame OXYZ in which the Sun and the Earth remain in fixed locations of the X axis, and we map the previous orbit to this reference frame. Finally, since the planetary perturbations introduce oscillations of amplitude d, we define the set ℬ as a neighborhood of the previous orbit of amplitude ρ ∈ (d, α/10), and we find orbits of the realistic model of the Solar System which remain in the set ℬ for a minimum of 200 years. For α < 10d, we apply the method using a ball ℬ centered at ℒ_{3} of radius α. For the minimum value of α = 2d we find an orbit stable up to 250 years in a neighborhood of ℒ_{3}. As an illustration of our method, we reproduce a sample of these orbits in Figs. 4 and 7.
Displacement X_{ℒ3} − X_{L3}, differences between the eigenvalues of ℒ_{3} and L_{3}, and amplitude of the dynamical substitute (DS) along the Y axis for each planet.
The paper is organized as follows: in Sect. 2 we introduce the MCR8BP and we study the system obtained by averaging the mean anomalies of the planets (except for the Earth); in particular, we introduce ℒ_{3} and compute the dynamical substitute in the original system conjugate to such equilibrium. In Sect. 3 we define the FLI method allowing us to recover orbits that remain in the neighborhood of L_{3} for at least 200 years in a model of the Solar System that is compatible with the JPL ephemerides.
2. The planar multicircular restricted planetary problem
2.1. Equations of motion
Let us consider the simplified planar model MCR8BP of our Solar System, where the eight planets from Mercury to Neptune perform circular motions around a fixed Sun in the same orbital plane. Our aim is to study the dynamics of a spacecraft gravitationally attracted by the planets close to the Sun–Earth Lagrangian equilibrium point L_{3}. We use μ and ε_{j} to denote the gravitational parameters of the Sun and of the jth planet, respectively, and n_{j} the mean motion of the jth planet. In an inertial reference frame Oxy centered at the Sun, the jth planet position is identified by r_{j}(t) = d_{j} (cos (n_{j}t + ϕ_{j}), sin (n_{j}t + ϕ_{j})), where d_{j} and ϕ_{j} denote the orbital radius and the initial value of the mean anomalies respectively. In order to compare the dynamics of the system with the dynamics of the Sun–Earth CR3BP we define the Hamiltonian of the system in a reference frame OXY centered at the Sun and such that the Earth position is (d_{3}, 0):
where
with j ≠ 3, where P_{X} and P_{Y} denote the momenta conjugate to X and Y, and , . We note that if ε_{j} = 0 for j ≠ 3, the Hamilton function (1) reduces to the Hamiltonian of the CR3BP; for this reason we treat ∑_{j ≠ 3}ε_{j}H_{j}(X, Y; t) as a quasiperiodic timedependent perturbation of the Sun–Earth CR3BP.
2.2. Averaging the planetary anomalies
In this section we introduce a timedependent near to the identity canonical transformation
conjugating the Hamiltonian (1) to a Hamiltonian of the form
where is of first order in the planetary masses, fourth order in the eccentricity, and of order where is defined as for j < 3 and for j > 3; ℛ″ is of second order in the planetary masses and
where is the elliptic function of first kind. In the following paragraph we show how the averaging Hamiltonian method allows us to compute the transformation 𝒞 explicitly. We note that and also ℛ″ are negligible with respect to for the time intervals of interest for astrodynamics. The dynamics of (2) are studied in Sect. 2.3.
The Hamiltonian averaging method. The starting point to average from the Hamiltonian (1) the dependence on the planetary anomalies (j ≠ 3) with a canonical transformation is to introduce the canonical modified Delaunay variables
as usual from the Delaunay variables L, G, l, g by
where
where a, e, M, and ω represent the osculating semimajor axis, eccentricity, mean anomaly, and argument of perihelion, and then we extend the phasespace by introducing the planetary mean anomalies θ_{j} and their conjugate momenta Θ_{j}. We obtain the Hamilton function
where
while 𝒲_{3} (Λ, Φ, λ, φ) and 𝒲_{j} (Λ, Φ, λ, φ, θ_{j}) are respectively H_{3} and H_{j} represented in the modified Delaunay variables. Because we are interested in the study of the dynamics close to L_{3}, we expand the perturbing function 𝒲_{j} in power series of the eccentricity up to order 3 and of or up to order 35 for j < 3 and j > 3, respectively. In this way, 𝒲_{j} can be represented as the Fourier series
As the semimajor axis associated to L_{3} is close to 1, when we represent the planetary perturbations we need to include terms up to to provide a sufficiently precise approximation: for example, with lower order approximations, the initial conditions provided by this theory would produce orbits that do not close to a periodic orbit like the ones represented in the top panels of Fig. 1.
Fig. 1.
Dynamical substitutes of ℒ_{3} in the Sun–Earth–Venus, Sun–Earth–Jupiter, Sun–Earth–Venus–Jupiter and complete systems obtained by integration. For the Sun–Earth–Venus–Jupiter and in the complete system the timespan is about 50 years. For all these cases we assume that for j ≠ 3. The red point indicates the position of L_{3}. 
In this context, it is natural to construct canonical transformations using the Lie method: given a generating function χ(Λ, Φ, λ, φ, θ), the Hamiltonian flow generated by χ at time t = 1 is a canonical transformation:
which can be represented by the Lie series:
where { , } denote the Poisson brackets, and conjugates the Hamiltonian 𝒲 to
The normal form 𝒲 does not depend on the planetary anomalies at first order in the planetary masses if we define the generating function by
with
where δ_{0} = 0, δ_{i} = 1 if i ≠ 0, and , , and . Precisely, the normal form corresponding to this generating function can be represented by:
with ε^{2} = ∑_{j ≠ 3}ε_{3}ε_{j} and
By introducing the averaged Cartesian coordinates ( ) from the averaged modified Delaunay variables (see Eq. (3)):
we have and
Finally, by representing and introducing , we have
2.3. Dynamics of the averaged Hamiltonian (2) close to ℒ_{3}
The Hamiltonian (2) admits an equilibrium point on the axis close to the Lagrangian point L_{3} of the CR3BP, which we denote by ℒ_{3}. More precisely, Hamiltonian (2) is a small perturbation of the Hamiltonian of the Sun–Earth CR3BP depending on . Therefore, on the one hand the equation:
is satisfied by any with ; on the other hand the equation
has a solution with X_{ℒ}_{3} close to X_{L3}. By including in the model all the planets from Mercury to Neptune we find X_{ℒ}_{3} − X_{ℒ}_{3} ∼ −2.58 10^{−7} AU ∼ −39 km. More precisely, by computing the individual contribution of each planet to the displacement of ℒ_{3} from L_{3} (see Table 1, where the individual contribution of a planet is computed by setting to zero the masses of all other planets, except for the Earth) we find that Venus and Jupiter each move the position of this equilibrium by about 218 and 176 km, respectively, in opposite directions; in the model where both the planets are included, their effects almost perfectly compensate for one another, leaving a displacement of about 40 km only. Table 1 shows that the individual contributions of the planets to the eigenvalues of the linearization matrix at ℒ_{3} are very small, and do not affect the nature of the equilibrium, which remains partially hyperbolic. As a consequence, around ℒ_{3} we find a family of periodic orbits, namely the Lyapunov orbits of ℒ_{3}.
2.4. The dynamical substitute: mapping ℒ_{3} back to the original Cartesian variables
The equilibrium point ℒ_{3} is finally mapped back to the original Cartesian variables P_{X}, P_{Y}, X, Y using the canonical transformation defined by the Hamiltonian flow at “time 1” of the generating function introduced in Eq. (10). Since the orbital elements of ℒ_{3} are characterized by a very small eccentricity (about 10^{−6}), it is mandatory to represent the generating function χ using the Poincaré variables , , , , where are related to by
Moreover, since for the times required by astrodynamics we are allowed to neglect in the transformation the terms proportional to the product of two planetary masses, we approximate the Lie series generated by χ with
Taking advantage of the small value of the eccentricity of ℒ_{3}, the numerical computation of the righthand sides of Eq. (11) can be done by neglecting the terms depending on at orders equal to or higher than four and for . Finally, we transform back to the Cartesian variables and we obtain the dynamical substitute of ℒ_{3} as a quasiperiodic function of time:
In the simpler bicircular fourbody problems, Sun–Earth–(jth planet), the dynamical substitute is a periodic function of time with frequency . In order to verify that there is a real orbit for the restricted problem that remains close to the dynamical substitute (X_{ℒ}_{3}(t),Y_{ℒ}_{3}(t)) for times that are of interest for astrodynamics, we numerically integrate the equations of motion of the restricted problem with initial conditions . In Fig. 1 we present the dynamical substitutes of ℒ_{3} computed by integrating the initial datum for four cases: the Sun–Earth–Venus, Sun–Earth–Jupiter, Sun–Earth–Venus–Jupiter, and the complete system (namely the system containing all eight planets) with for all j ≠ 3. In Table 1 we also present the libration amplitudes for each restricted Sun–Earth–(jth planet) problem. We find that in the MCR8BP, the orbit is a quasiperiodic oscillation around L_{3} with amplitude d ∼ 10^{−4} AU. Each orbit in a neighborhood of ℒ_{3} can be mapped back to the original Cartesian variables using the canonical transformations (11) and 𝒟; for example, the periodic Lyapunov orbits will be mapped back to quasiperiodic orbits of amplitudes larger than d. Therefore, in the MCR8BP we have a family of quasiperiodic oscillations around L_{3} of amplitudes α > d.
3. The fast Lyapunov indicator method and the computation of the librations in a realistic model of the Solar System
As the restricted multicircular model represents just a simplified model of our Solar System, we must investigate whether or not orbits close to L_{3} still exist when the trajectories of the planets are realistic. Furthermore, as the effects of Mercury on both the dynamical substitute and in the change of the other planets ephemerides are negligible for the short time intervals that we are taking into account, we consider a model of the Solar System that includes the planets from Venus to Neptune. We use r = (x, y, z) and v = (v_{x}, v_{y}, v_{z}) to denote the position and velocity vectors of the spacecraft, and r_{j} = (x_{j}, y_{j}, z_{j}) and v_{j} = (v_{x, j}, v_{y, j}, v_{z, j}) the position and velocity vector of the jth planet in the Heliocentric reference frame. To identify a configuration generalizing the Lagrangian point L_{3} of the CR3BP in this more realistic system, we first introduce a Heliocentric rotatingpulsating reference frame OXYZ, whose Cartesian axes are parallel, at any time t, to the unit vectors
In a such reference frame, the Earth state vector is (1, 0, 0, 0, 0, 0) at any time t. We apply the FLI method in order to find orbits staying close to L_{3} in the realistic model for centuries; this method provides initial conditions of trajectories remaining in a given neighborhood ℬ of a libration orbit found in an approximated model. For the definition of ℬ we have to consider both the amplitude of the libration orbit and the differences between the realistic model and the CR3BP due to Earth’s elliptic orbit. Hence, for the application of the FLI method we first compute a planar Lyapunov orbit of the CR3BP with amplitude α and then we use this orbit to find librations in the ER3BP remaining close to it. We denote these librations by ℓ(α) and the details of their computation are explained in Sect. 3.2. By assuming that the planetary perturbations introduce additional oscillations of amplitude d comparable to those found in the MCR8BP, we finally look for hyperbolic orbits of the realistic model which remain in a neighborhood of radius 2d of the reference orbit ℓ(α) for centuries (Sect. 3.3). We compute both ℓ(α) of the ER3BP for different orders of the amplitude and the libration orbits in the realistic model of the Solar System using the recent version of the modified FLI method.
3.1. The FLI method
The FLI method was originally introduced by Froeschlé et al. (1997a,b) to numerically distinguish between regular and chaotic orbits. In recent years, the method has been modified in order to compute hyperbolic orbits of dynamical systems and their asymptotic solutions: Guzzo & Lega (2014, 2018), and Lega & Guzzo (2016) provide a detailed theoretical justification of the method, with a focus on the computation of transit orbits of the SunJupiter system in the circular restricted threebody problem; in Guzzo & Lega (2015, 2017) the method has been adapted to compute the collision manifold of the comet 67P/ChuryumovGerasimenko in a model of the Solar System that is compatible with the JPL ephemerides. In this paper, for the first time, the modified version of the FLI method is used to investigate the dynamics originating at a Lagrangian equilibrium point in a realistic model of the Solar System. Moreover, since we are considering the point L_{3} of the Sun–Earth system, the CR3BP provides only a very crude approximation. The extension from the CR3BP problem to the realistic one required a new strategy based on the preliminary study of intermediate models to include the planetary perturbations and the eccentricity of Earth’s orbit. In this way we first find the orbits ℓ(α) of the Sun–Earth ER3BP, and then the libration orbits of the realistic model of the Solar System. In particular, this usage of the chaos indicator provides an innovative and novel application of dynamical systems theory to astrodynamics, where the shortperiod perturbations represent a relevant part of the model.
The modified FLI is introduced as follows. We first consider the spacecraft equations of motion written as a system of firstorder differential equations:
where ξ denotes the state vector of the spacecraft in a given reference frame. For simplicity, we choose as epoch t = 0 the time of an Earth passage to the perihelion. The variational equations of (13) are
where denotes the Jacobian matrix of F. If we aim to find the hyperbolic orbits that remain in a neighborhood ℬ of a given target orbit γ for as long as possible, we introduce the window function w(R) and the modified FLI defined as
respectively, where R(s) is the position vector of the spacecraft in the rotatingpulsating reference frame computed by integrating the differential equations (13) with initial conditions ξ_{0}, and Ξ(s) is the solution of the variational equations (14) with Ξ_{0} as initial tangent vector. When the target orbit γ is partially hyperbolic and ℬ is suitably small, the only possibility for an orbit to increase the modified FLI to the highest possible values is to remain very close to the target over the whole time interval [0, T]. Here, we have the two possible cases: (i) a large measure set ℬ_{I} ⊂ ℬ of initial conditions corresponds to solutions that exit from ℬ in short times T_{0} and whose mFLI does not increase for T > T_{0}; (ii) a small measure set ℬ_{II} ⊂ ℬ corresponds to solutions which do not exit from ℬ, remain close to a partially hyperbolic orbit of ℬ, and whose mFLI increases linearly with time T (Guzzo & Lega 2014). Orbits that come back to (or enter) the set ℬ, but do not stay close to the target over the time interval of interest [0, T] necessarily enter the set ℬ_{I} and exit from ℬ in a short time T_{0}. Therefore, the modified FLI sharply identifies the orbits (II) when it is computed on grids of initial conditions that are transverse to ℬ_{II}. The key factor for the success of the method is the suitable choice of the set ℬ and of the grid of initial conditions, which is suggested by the application of perturbation theories and approximated models.
3.2. Computation of libration orbits in the ER3BP
The Sun–Earth ER3BP deals with the study of the spacecraft dynamics gravitationally attracted by the Sun and the Earth while the latter performs an elliptic motion around the barycenter of the Sun–Earth system; in order to compare the ER3BP with the CR3BP, the equations of motion of the ER3BP are typically described in the rotatingpulsating reference frame (Szebehely 1967). In order to find the orbits ℓ(α) of the ER3BP, we use the modified FLI (15) using a planar Lyapunov orbit of the CR3BP of libration amplitude α as target orbit γ. The set ℬ is a neighborhood of γ with amplitude ; if (X_{0}, 0, 0, 0, V_{Y, 0}, 0) is a point of the target Lyapunov orbit, we choose the grid of initial data as the set (X_{0}, 0, 0, 0, V_{Y}, 0), where V_{Y} is evenly spaced over some interval. The extrema of this latter interval are chosen such that it contains a strict maximum of the modified FLI. We then restrict this interval around the maximum and repeat the computation of the FLI for the new grid and for a longer integration time. The procedure is iterated until we achieve a sufficiently long integration time. In Fig. 2 we show the results of the method for a timespan of 1000 years and for values of α equal to 3 × 10^{−2}, 10^{−3}, and 10^{−4} AU, correspondingly to equal to 10^{−3}, 7 × 10^{−5}, and 10^{−5} AU. We note that such orbits are slowly drifting along the Y axis (at a speed of about 2 × 10^{−6}, 5 × 10^{−8} and 4.5 × 10^{−9} AU yr^{−1} for ℓ(3 × 10^{−2} AU), ℓ(10^{−3} AU), and ℓ(10^{−4} AU), respectively): these drifts can be reduced by further refining the interval of the initial data.
Fig. 2.
Representation of ℓ(α) orbits (gray lines) of the ER3BP in the rotatingpulsating reference frame with α = 3 × 10^{−2}, 10^{−3}, 10^{−4} AU for 1000 years. The values of in the three different cases are 10^{−3}, 7 × 10^{−5}, and 10^{−5} AU. The representation of ℓ(α) in the timespan of 1 year (black curves) emphasizes the small drifts of the librations along the Y axis. The red point indicates the position of L_{3}. 
3.3. Computation of the libration orbits in the realistic model of the Solar System
In order to study the spacecraft dynamics in a model whose ephemerides are comparable with those provided by the JPL computation service, we consider the spacecraft equations of motion in the barycentric reference frame in a model composed of the Sun and seven planets whose initial conditions are provided by JPL (for the Earth we consider the Earth–Moon barycenter). To illustrate how to find libration orbits in such a realistic model, we apply the FLI method and we consider ℓ(3 × 10^{−2} AU) and ℓ(10^{−3} AU) as target orbits and take ℬ as a neighborhood of ℓ(3 × 10^{−2} AU) and ℓ(10^{−3} AU) of amplitude equal to 10^{−3} and 2 × 10^{−4} AU, respectively. The computation of the modified FLI is stopped once we have restricted the set of initial data such that the interval along the V_{Y} is smaller than 5 × 10^{−9} AU day^{−1} and once we achieve an integration time of centuries (in the cases of the librations with amplitudes of 3 × 10^{−2} and 10^{−3} AU the timespans for the modified FLI computation are 1000 and 400 years respectively). When we conclude this iterative method, we obtain the modified FLI profile, in which the maximum is sharply identified by the initial datum . For example, the modified FLI profile for the research of the libration with amplitude 3 × 10^{−2} AU is represented in Fig. 3. The trajectory associated to the initial datum with the maximum mFLI value (blue point in Fig. 3) represents the orbit that remains the closest to the target one. Moreover, this trajectory (see the first two panels in Fig. 4) is librating around L_{3} with an amplitude of 3 × 10^{−2} AU. To appreciate the precision of the method, we also considered the orbits associated to another two initial data with values of the modified FLI that are different from the maximum (the magenta and brown points in Fig. 3). Such orbits, shown in Fig. 5, are not librations around L_{3}, even if the variation of the initial velocity is very small (about 3 × 10^{−12} and −3 × 10^{−11} AU day^{−1}, respectively). Now, we focus our attention on the orbits associated to the maximum value of mFLI for both the case of ℓ(3 × 10^{−2} AU) and ℓ(10^{−3} AU) as target orbits. These trajectories are shown in Fig. 4 in the timespans of 650 and 200 years respectively. In the same figure we also show the representations of such trajectories restricted to a timespan of 20 years (black curves) which highlight the drift along the Y axis (about 4 × 10^{−6} AU yr^{−1}) caused by the planetary perturbations. The relative drift is greater for the libration with amplitude 10^{−3} AU and is smaller for the orbit with greater amplitude. Hence, the orbit with the greater amplitude remains closer to its target orbit. Moreover, the value of for the orbit with smaller amplitude is close to the amplitude of the oscillation introduced by the planetary perturbations (last panel in Fig. 1): this implies that values of are too small and for such a choice of the FLI method does not converge to a libration orbit around L_{3}. Therefore, to find oscillations around L_{3} with amplitudes smaller than 2d it is necessary to apply the FLI method by using a sphere in (X, Y, Z) centered at L_{3} with radius as ℬ. We note also that the application of the modified FLI method with a onedimensional grid of initial data is not sufficient for the research of librations because of the degrees of freedom of the problem: we have to use a twodimensional grid of initial data (X, 0, 0, 0, V_{Y}, 0) centered at L_{3}. Therefore, we applied this method with and AU and found (see e.g., Fig. 6) that the bigger values of the mFLI lie on a line not parallel to the X and V_{Y} axis. This complicates the procedure of refining the initial data, which is summarized in the caption of the figure. Finally, Fig. 7 shows the librations obtained using the twodimensional grid with and AU: the timespans of these two orbits are 700 and 250 years respectively. Hence, the FLI method allows us to find libration orbits of the realistic models with amplitudes bigger than d (up to 0.03 AU) for timespans of centuries.
Fig. 3.
Profile of the modified FLI for the timespan of 1000 years in the case of the research of the libration with amplitude 3 × 10^{−2} AU in the realistic model. The blue point denotes the location of , namely the velocity component along Y characterized by the maximum mFLI value, and the trajectory associated to the initial condition is shown in the first two panels of Fig. 4. We also show the locations of another two initial conditions (the magenta and brown point) with a different value of the chaos indicator whose trajectories are plotted in Fig. 5. 
Fig. 4.
Representation of the librations around L_{3} in the realistic model obtained by changing the target orbit ℓ(α) and . The first two figures represent the librations found by using ℓ(3 × 10^{−2} AU) and ℓ(10^{−3} AU), respectively. The black lines represent the restriction of the trajectories to 20 years. The red point indicates the position of L_{3} in the rotatingpulsating reference frame. 
Fig. 5.
Trajectories associated to the magenta and brown point in Fig. 3. The timespan for both orbits is 650 years. The plots highlight the precision of the modified FLI method. Small displacements of the initial conditions corresponding to the maximum mFLI value lead to large differences in the orbits: in this scenario, the variation of the initial velocity is about 3 × 10^{−12} and −3 × 10^{−11} AU day^{−1} for the magenta and brown trajectories. Both orbits are escaping from L_{3}. 
Fig. 6.
Representation of the modified FLI map associated to the bidimensional grid of initial data in the case in which L_{3} is the reference point and ; the colors indicate the values of the chaos indicator computed in the timespan of 1000 years. We note that the maximum values of mFLI are not in a line parallel to the X or V_{Y} axis. For this reason, if we want to have a good refinement of the twodimensional grid, we need to isolate two points P_{0} = (X_{0}, V_{Y, 0}), P_{1} = (X_{1}, V_{Y, 1}) (green points) of the twodimensional grid characterized by a high value of the modified FLI and then take points X and V_{Y} such that X = X_{0} + δ(X_{1} − X_{0})+ϵ(V_{Y, 1} − V_{Y, 0}), V_{Y} = V_{Y, 0} + δ(V_{Y, 1} − V_{Y, 0})−ϵ(X_{1} − X_{0}) with δ ∈ [0, 1] and ϵ ∈ [ − Δϵ, Δϵ]: we used Δϵ equal to 2 × 10^{−5} for the refinement of the grid shown in this figure. 
Fig. 7.
Representation of the librations around L_{3} in the realistic model obtained assuming a neighborhood of L_{3} with amplitude equal to 10^{−3} and 2 × 10^{−4} AU, respectively. The timespans of the two orbits are 700 and 250 years respectively. The black lines represent the restriction of the trajectories to 20 years. The red point indicates the position of L_{3} in the rotatingpulsating reference frame. 
4. Conclusions
Since close to the Sun–Earth Lagrangian point L_{3} the planetary perturbations (mainly from Jupiter and Venus) are stronger than Earth’s attraction, it is not evident that the dynamics close to L_{3} defined within the circular restricted threebody problem can be continued to a realistic model of the Solar System. We analyze in particular the effects of the planetary perturbations and of the eccentricity of the Earth. The planetary perturbations provide oscillations of about 15 000 km, while the effect of the eccentricity of the Earth, apart from a small oscillation, is absorbed by considering a rotatingpulsating frame. Using a combination of the Hamiltonian averaging method with the FLI method we find orbits that remain close to L_{3} for at least 200 years in a model of the Solar System whose dynamics are compatible with JPL digital ephemerides, with libration amplitudes ranging from 2.5 × 10^{−4} to 3 × 10^{−2} AU.
Acknowledgments
We thank the anonymous referee for useful remarks and suggestions. M.G. acknowledges the project MIURPRIN 20178CJA2B “New frontiers of Celestial Mechanics: theory and applications”.
References
 Barrabés, E., & Ollè, M. 2006, Nonlinearity, 19, 2065 [CrossRef] [Google Scholar]
 Froeschlé, C., Gonczi, R., & Lega, E. 1997a, Planet. Space Sci., 45, 881 [NASA ADS] [CrossRef] [Google Scholar]
 Froeschlé, C., Lega, E., & Gonczi, R. 1997b, Celestial Mech. Dyn. Astron., 67, 41 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Gómez, G., & Mondelo, J. 2001, Physica D, 157, 283 [NASA ADS] [CrossRef] [Google Scholar]
 Gómez, G., Llibre, J., Martínez, R., & Simó, C. 2001, in Dynamics and Mission Design Near Libration Points, World Scientific Monograph Series in Mathematics (World Scientific Publishing Co. Inc.), 1 [Google Scholar]
 Gómez, G., Koon, W., Lo, M., et al. 2004, Nonlinearity, 17, 1571 [NASA ADS] [CrossRef] [Google Scholar]
 Guzzo, M., & Lega, E. 2013, MNRAS, 428, 2688 [NASA ADS] [CrossRef] [Google Scholar]
 Guzzo, M., & Lega, E. 2014, Soc. Ind. Appl. Math., 74, 1058 [CrossRef] [Google Scholar]
 Guzzo, M., & Lega, E. 2015, A&A, 579, 1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Guzzo, M., & Lega, E. 2017, MNRAS, 469, S321 [CrossRef] [Google Scholar]
 Guzzo, M., & Lega, E. 2018, Physica D, 373, 38 [CrossRef] [Google Scholar]
 Hou, X., Tang, J., & Liu, L. 2007, NASA Technical Report: 20080012700 [Google Scholar]
 Jorba, À., & Nicolás, B. 2020, Communications in Nonlinear Science and Numerical Simulation, 89, 105327 [CrossRef] [Google Scholar]
 Koon, W. S., Lo, M. W., Marsden, J. E., & Ross, S. D. 2006, Dynamical Systems, The ThreeBody Problem, and Space Mission Design (Pasadena, CA, USA: California Institute of Technology) [Google Scholar]
 Lega, E., & Guzzo, M. 2016, Physica D, 325, 41 [CrossRef] [Google Scholar]
 Lega, E., Guzzo, M., & Froeschlé, C. 2011, MNRAS, 418, 107 [NASA ADS] [CrossRef] [Google Scholar]
 Martin, C., Conway, B. A., & Ibánez, P. 2010, in Space Manifold Dynamics: Novel Spaceways for Science and Exploration, eds. E. Perozzi, & S. FerrazMello (Springer) [Google Scholar]
 Pàez, R., & Efthymiopoulos, C. 2015, Celestial Mech. Dyn. Astron., 121, 139 [Google Scholar]
 Parker, J., & Anderson, R. L. 2013, LowEnergy Lunar Trajectory Design (Pasadena, CA: Jet Propulsion Laboratory, California Institute of Technology) [Google Scholar]
 Szebehely, V. 1967, Theory of Orbits: The Restricted Problem of Three Bodies (Academic Press) [Google Scholar]
 Tantardini, M., Fantino, E., Ren, Y., et al. 2010, Celestial Mech. Dyn. Astron., 108, 215 [CrossRef] [Google Scholar]
 Terra, M. O., Simó, C., & de Sousa Silva, P. 2014, Proc. IAC, 6, 4535 [Google Scholar]
 van Damme, C. C., Gorgojo, R. C., GilFernandez, J., & Graziano, M. 2010, in Space Manifold Dynamics: Novel Spaceways for Science and Exploration, eds. E. Perozzi, & S. FerrazMello (Springer) [Google Scholar]
All Tables
Displacement X_{ℒ3} − X_{L3}, differences between the eigenvalues of ℒ_{3} and L_{3}, and amplitude of the dynamical substitute (DS) along the Y axis for each planet.
All Figures
Fig. 1.
Dynamical substitutes of ℒ_{3} in the Sun–Earth–Venus, Sun–Earth–Jupiter, Sun–Earth–Venus–Jupiter and complete systems obtained by integration. For the Sun–Earth–Venus–Jupiter and in the complete system the timespan is about 50 years. For all these cases we assume that for j ≠ 3. The red point indicates the position of L_{3}. 

In the text 
Fig. 2.
Representation of ℓ(α) orbits (gray lines) of the ER3BP in the rotatingpulsating reference frame with α = 3 × 10^{−2}, 10^{−3}, 10^{−4} AU for 1000 years. The values of in the three different cases are 10^{−3}, 7 × 10^{−5}, and 10^{−5} AU. The representation of ℓ(α) in the timespan of 1 year (black curves) emphasizes the small drifts of the librations along the Y axis. The red point indicates the position of L_{3}. 

In the text 
Fig. 3.
Profile of the modified FLI for the timespan of 1000 years in the case of the research of the libration with amplitude 3 × 10^{−2} AU in the realistic model. The blue point denotes the location of , namely the velocity component along Y characterized by the maximum mFLI value, and the trajectory associated to the initial condition is shown in the first two panels of Fig. 4. We also show the locations of another two initial conditions (the magenta and brown point) with a different value of the chaos indicator whose trajectories are plotted in Fig. 5. 

In the text 
Fig. 4.
Representation of the librations around L_{3} in the realistic model obtained by changing the target orbit ℓ(α) and . The first two figures represent the librations found by using ℓ(3 × 10^{−2} AU) and ℓ(10^{−3} AU), respectively. The black lines represent the restriction of the trajectories to 20 years. The red point indicates the position of L_{3} in the rotatingpulsating reference frame. 

In the text 
Fig. 5.
Trajectories associated to the magenta and brown point in Fig. 3. The timespan for both orbits is 650 years. The plots highlight the precision of the modified FLI method. Small displacements of the initial conditions corresponding to the maximum mFLI value lead to large differences in the orbits: in this scenario, the variation of the initial velocity is about 3 × 10^{−12} and −3 × 10^{−11} AU day^{−1} for the magenta and brown trajectories. Both orbits are escaping from L_{3}. 

In the text 
Fig. 6.
Representation of the modified FLI map associated to the bidimensional grid of initial data in the case in which L_{3} is the reference point and ; the colors indicate the values of the chaos indicator computed in the timespan of 1000 years. We note that the maximum values of mFLI are not in a line parallel to the X or V_{Y} axis. For this reason, if we want to have a good refinement of the twodimensional grid, we need to isolate two points P_{0} = (X_{0}, V_{Y, 0}), P_{1} = (X_{1}, V_{Y, 1}) (green points) of the twodimensional grid characterized by a high value of the modified FLI and then take points X and V_{Y} such that X = X_{0} + δ(X_{1} − X_{0})+ϵ(V_{Y, 1} − V_{Y, 0}), V_{Y} = V_{Y, 0} + δ(V_{Y, 1} − V_{Y, 0})−ϵ(X_{1} − X_{0}) with δ ∈ [0, 1] and ϵ ∈ [ − Δϵ, Δϵ]: we used Δϵ equal to 2 × 10^{−5} for the refinement of the grid shown in this figure. 

In the text 
Fig. 7.
Representation of the librations around L_{3} in the realistic model obtained assuming a neighborhood of L_{3} with amplitude equal to 10^{−3} and 2 × 10^{−4} AU, respectively. The timespans of the two orbits are 700 and 250 years respectively. The black lines represent the restriction of the trajectories to 20 years. The red point indicates the position of L_{3} in the rotatingpulsating reference frame. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.