Revisiting the averaged problem in the case of mean-motion resonances in the restricted three-body problem

A classical approach to the restricted three-body problem is to analyze the dynamics of the massless body in the synodic reference frame. A different approach is represented by the perturbative treatment: in particular the averaged problem of a mean-motion resonance allows to investigate the long-term behavior of the solutions through a suitable approximation that focuses on a particular region of the phase space. In this paper, we intend to bridge a gap between the two approaches in the specific case of mean-motion resonant dynamics, establish the limit of validity of the averaged problem and take advantage of its results in order to compute trajectories in the synodic reference frame. After the description of each approach, we develop a rigorous treatment of the averaging process, estimate the size of the transformation and prove that the averaged problem is a suitable approximation of the restricted three-body problem as long as the solutions are located outside the Hill’s sphere of the secondary. In such a case, a rigorous theorem of stability over finite but large timescales can be proven. We establish that a solution of the averaged problem provides an accurate approximation of the trajectories on the synodic reference frame within a finite time that depend on the minimal distance to the Hill’s sphere of the secondary. The last part of this work is devoted to the co-orbital motion (i.e., the dynamics in 1:1 mean-motion resonance) in the circular-planar case. In this case, an interpretation of the solutions of the averaged problem in the synodic reference frame is detailed and a method that allows to compute co-orbital trajectories is displayed.


Introduction
This work focuses on the restricted three-body problem, that is the study of the motion of a massless body affected by the gravitational attraction of two massive bodies. More precisely, we will consider the situation for which the mass of the secondary body is treated as a small quantity. Since the planetary three-body problem will also be mentioned, we recall that it corresponds to the study of the motion of two massive bodies orbiting a more massive one, the three bodies being governed only by their mutual gravitational interactions.
The analysis of the dynamics in the synodic reference frame, that is the frame rotating with the mean longitude of the secondary, is the classical approach adopted for the restricted three-body problem. Usually, periodic orbit families and the dynamics located in their neighborhood are computed by using Poincaré maps and continuation methods (see, e.g., [13,37]).
Perturbative treatments provide another approach. They allow to investigate specific regions of the phase space through a proper approximation. Among them, averaging methods are common techniques in order to study the long-term dynamics of the solutions. For instance, the secular problem studies the long-term deformation of the ellipse of the massless body as well as the evolution of its orientation in the tridimensional space. It is obtained by the averaging of the Hamiltonian over the mean longitudes of the secondary and of the massless body. More precisely, it corresponds to a symplectic transformation, that is supposed to be close to the identity, and that maps the original Hamiltonian to the secular one.
Lagrange [17] introduced the secular problem in the framework of the stability of the Solar System and the expression of the secular Hamiltonian of the planetary three-body problem was given by Poincaré [26]. Precise estimates on the size of the transformation of averaging were required in order to prove theorems of stability like KAM theory and were provided, especially by Arnol'd [1], Féjoz [10], Chierchia and Pinzari [6].
When the massless body is in mean-motion resonance with the secondary, that is, when their orbital periods are commensurable, the transformation leading to the secular Hamiltonian is no more close to the identity and the solutions of the secular problem do not provide a good representation of the real motion. In such a case, it is still possible to use averaging techniques: the averaging process is performed over one mean longitude, generally the one of the secondary, and after the introduction of a resonant angle, that is a particular linear combination of the two mean longitudes which characterizes the mean-motion resonance. This defines the averaged problem that will be considered in this work.
Many authors investigated mean-motion resonances through an averaged Hamiltonian and the literature on this subject has become so rich that it is impossible to cite all the articles here. Nevertheless, let us mention the important series of works realized by Schubart [31][32][33], which took advantage of the canonical variables and method suggested by Poincaré [27] and applied an averaging process in order to get the interesting part of the Hamiltonian for mean-motion resonances. Likewise, the second fundamental model of resonance, developed by Henrard and Lemaitre [16], follows the strategy of Poincaré [27], and is commonly used in order to study mean-motion resonances. Moons [19] extended the work of Schubart and presented an integrator adapted to the solutions of the averaged problem. Being the latter not valid for the 1:1 mean-motion resonance, Nesvorný et al. [24] adapted the algorithm with a different choice of canonical variables.
The co-orbital motion, or equivalently, the trajectories in 1:1 mean-motion resonance with the secondary, has been intensively studied in the framework of the averaged problem (see, e.g., [18,20,34,35]). In such a case, since the semimajor axis of the massless body is almost the same as the one of the secondary, the issue generated by periodical close encounters arises, even for quasi-circular trajectories. In particular, Robutel and Pousse [30] and Pousse et al. [28] highlighted, with the help of a frequency analysis, that the averaged Hamiltonian reflects poorly the dynamics close to the singularity associated with the collision between the secondary and the massless body. Rigorous estimates on the averaging process have been given by Robutel et al. [29]. More precisely, they allowed the authors to prove in the planetary three-body problem, that the averaged problem is valid for two co-orbital bodies on quasi-circular orbits that stand at a mutual distance larger than their respective Hill's radius.
The limit of validity of the averaged problem is not specific to the case of the co-orbital motion and can also occur for other resonant trajectories that cross the orbit of the secondary (i.e., trajectories with a non-negligible eccentricity). This weakness was already outlined in the works of Schubart [31] and Moons [19]. Therefore, in the present paper, we intend to generalize the result given by Robutel et al. [29] and provide rigorous estimates on the averaging process in order to define a domain of validity of the averaged problem, in the case of a generic mean-motion resonance, and for any value of inclination and eccentricities (massless body as well as secondary).
According to the Poincaré classification (see, e.g., [5] for more details), some of the periodic families described in the synodic reference frame are related to mean-motion resonances and thus can also be tackled in the averaged problem as defined here. For that reason, we also intend to bridge a gap between the classical approach in the synodic reference frame and the aver-aged problem with a unified Hamiltonian formalism that allows to represent solutions in both approaches. The underlying idea of this work is to understand the limit of validity of the averaged problem and take advantage of its solutions (e.g., initial conditions, types of motion, frequencies) for the computation of trajectories in the synodic reference frame.
The paper is structured as follows. Section 2 introduces the restricted three-body problem through the classical approach, recalls some remarkable solutions in the synodic reference frame that are connected to the co-orbital motion, and presents the reasoning that led to the averaged problem.
In Sect. 3, the size of the transformation of averaging is estimated. This allows to define a domain of validity of the approximation and to prove a rigorous theorem of stability over finite times.
In Sect. 4, we focus on the co-orbital motion in the circular-planar case and detail the correspondence between a solution of the averaged problem and its corresponding trajectory in the synodic reference frame. In particular, we will recover the remarkable solutions described in Sect. 2. Finally, a method that allows to compute co-orbital trajectories in the synodic reference frame will be described.
Appendix A gives the proof of the theorems and lemma used in our reasonings.

Definition in the heliocentric reference frame
Let (r,ṙ) be, respectively, the heliocentric position and velocity vector in R 3 of a massless body (particle, spacecraft or asteroid), that is affected by the gravitational attraction of a massive primary (the Sun or a planet) of mass 1 − ε > 1/2, and a secondary (a planet or a moon) of mass ε > 0. The motion of the two massive bodies, respectively, denoted as Sun and planet, follows a solution of the two-body problem. Hence, the trajectory of the planet, denoted r (t) in the heliocentric reference frame, lies on an ellipse that can be defined by the orbital elements (a , e , I , Ω , ω , v ), i.e., respectively, semimajor axis, eccentricity, inclination, longitude of the node, argument of the periaster and true anomaly. Without loss of generality, the scale and orientation of the orbit are arbitrarily chosen such that (a , I , Ω , ω ) = (1, 0, 0, 0).
Likewise, the orbital period of the planet is fixed to 2π (and therefore its mean motion to 1) which imposes the gravitational constant to be equal to 1. The eccentricity e is a parameter of the problem associated with the shape of the planet's orbit while the angle v stands for its position on the ellipse. Instead of using the true anomaly, the mean longitude λ will be adopted in order to take advantage of its proportionality to time since In the heliocentric reference frame, the equations of motion of the particle read (see [22]): where " · " is the Euclidean norm associated with the scalar product denoted " • " in what follows. The two first terms are, respectively, the gravitational force of the Sun and of the planet. The third term is associated with the acceleration of the heliocentric reference frame generated by the Sun-planet gravitational interactions.

Hamiltonian formalism
Since the heliocentric vectors r andṙ are canonical variables, then the Hamiltonian function provides the equations of motion (1). As H depends on the periodicity of the planet, it is non-autonomous. Moreover, the system just written, that describes the dynamics of the particle, has 3 degrees of freedom associated with the position and velocity vectors in the tridimensional space. A classical technique that allows to overcome the non-autonomous character of the system consists in extending the phase space with the addition of a generic variableΞ conjugated with λ . In this extended phase space, the Hamiltonian reads H +Ξ , it has 4 degrees of freedom and describes the coupled motion of the particle and the planet, namely As a consequence, investigating the restricted threebody problem consists in studying an autonomous ODE whose solutions belong to a 8-dimensional phase space (position and velocity in the tridimensional space, the mean longitude of the planet and its conjugated variable). A classical approach in order to simplify the analysis is to consider the behavior of the particle in the synodic reference frame, that is the frame that rotates with λ in the orbital plane of the planet.
2.2 A classical approach: the synodic reference frame

The synodic reference frame
Let us denote R k (α), the rotation matrix of an angle α about the k-axis (k ∈ {1, 2, 3}), and L(r,ṙ) = r ×ṙ, the angular momentum of the particle in the heliocentric reference frame. We recall that due to the influence of the planet, L(r,ṙ) is not a conserved quantity of the restricted three-body problem.
With the help of the Hamiltonian formalism, the symplectic transformation associated with the synodic reference frame reads

It provides the Hamiltonian
with R (λ ) = R 3 (−λ )r (λ ) that corresponds to the position of the planet. Moreover, the velocity of the particle in the synodic reference frame is deduced as follows: This framework gives rise to an important reduction in the context of the circular case (e = 0): since the planet is a fixed point located in R = (1, 0, 0), the Hamiltonian does not depend on λ and Ξ SF is an integral of motion that can be dropped. Hence, the dimension of the phase space to explore is reduced by two units. Moreover, H SF is related to the Jacobi constant that defines the energy of the particle. In the synodic reference frame centered on the Sun, the Jacobi constant can be written as follows: Thus, for a given value of C , the corresponding isoenergetic hypersurface is a manifold of dimension 5. C being the only conserved quantity (see [37]), it is not possible to reduce the problem through another global transformation. A further way to simplify the study is to consider the particle's motion restricted to the orbital plane of the planet. Hence, the phase space to explore can be reduced by two units. Consequently, investigating the restricted three-body problem in the circularplanar case is equivalent to explore a one-parameter family of 3-dimensional manifolds parametrized by the energy. Without too much details, the following part is dedicated to some remarkable solutions of the circularplanar case that are relevant for the scope of this work, i.e., the families that are connected to the co-orbital motion.

Some remarkable solutions in the circular-planar case (e = 0)
We recall that the configuration space of the circularplanar case coincides with the orbital plane of the planet. In the following, the motion of the particle will be described in terms of polar coordinates with φ = arg(R) that illustrates the relative motion between the planet and the particle and R = R . Most of the results mentioned below can be found with different notations in the book of Szebehely [37].  1 Periodic orbits in the synodic reference frame for a Sun-Jupiter-like system (ε = 1/1000) in the circular-planar case. More precisely, the periodic orbits belong to (a) the Lyapunov family L 3 , (b) the short-periodic family L s 4 , (c) the long-periodic family L l 4 , and (d) the family f . Their initial conditions are computed with the help of the Poincaré maps defined by the following sections: First of all, the five Lagrange fixed points, denoted L j for j = 1, 2, 3, 4, 5, are the unique equilibria of the restricted three-body problem in the circular case. L 1 and L 2 belong to the Sun-planet axis, in that is, from either side of the planet. Moreover, they embody the diameter of the Hill's sphere of the planet, that is the region of the configuration space inside which the gravitational influence of the planet dominates with respect to the one of the Sun. L 3 is also located on the Sun-planet axis, in L 4 and L 5 are the Lagrange configurations such that the particle lie at the vertex of an equilateral triangle formed with the Sun and the planet, that is, in φ j = (−1) j × 60 • and R j = 1.
For j = 4, 5 and ε small enough 1 , L j is an elliptic equilibrium where two one-parameter families of periodic orbits stem from. They are tangential to each center eigenspace of the equilibrium point. Being the two center eigenspaces associated with frequencies, respectively, in O(1) and O( √ ε), these families are generally denoted as short-periodic L s j and long-periodic L l j , in correspondence to their associated timescale in the neighborhood of the equilibrium. L 1 , L 2 and L 3 are unstable for all ε > 0 and each equilibrium possesses one center eigenspace. The same reasoning applies and provides three one-parameter families of periodic orbits generally known as the Lyapunov families L 1 , L 2 and L 3 . Only L 3 will be discussed in the following.
The Poincaré map is the classical way to compute periodic orbits. For Considering a periodic orbit that belongs to each family (black curve) whose cross section is denoted (X 0 , Y 0 ,Ẋ 0 ,Ẏ 0 ) a trajectory located in its vicinity, that is, with an initial condition (X, Y,Ẋ ,Ẏ ) = tions are given by Σ = {Y = R j sin φ j ,Ẏ > 0} and Σ ∩ {X > 0}, which require three free parameters (e.g., the energy, X andẊ ) in order to locate the crossing. We recall that the Lyapunov trajectories are symmetrical with respect to the Sun-planet axis and cross the X -axis inẊ = 0. Thus, a natural section, that requires only two parameters (e.g., the energy and X ), is given by Σ 0 = {Y = 0,Ẏ < 0,Ẋ = 0}. Then, a fixed-point method is generally performed from a suitable initial guess that makes the method convergent. For that purpose, a first approximation of a crossing is obtained by the resolution of the linearized system associated with the equilibrium L j and a continuation method is implemented. Figure 1 displays some periodic orbits, computed in the case of a Sun-Jupiter-like system (ε = 1/1000) by varying X along the section.
Since X > −R 3 increases, the size of a trajectory that belongs to L 3 (Fig. 1a) increases and its shape no longer looks like to an ellipse centered on L 3 . More precisely, φ and R oscillate, respectively, about 180 • and 1, whose respective amplitude increases with X and reaches large values close to 180 • and 1. Moreover, the "guiding center" of each periodic trajectory (the approximate position around which the trajectory oscillates) remains L 3 . Decreasing X < 1/2, the shape of the trajectories of L s 4 (Fig. 1b) has a quite similar evolution to the one of L 3 . However, two main differences exist: the shape is not symmetric, and the guiding center shifts from L 4 toward L 3 along the circle R = 1. The same behavior is observed symmetrically for L s 5 . We point out that, for a given energy, L s 4 and L s 5 merge together with L 3 . This result was found by Deprit et al. [9] for an Earth-Moon-like system (ε = 1/81) in the circular-planar case. The features of L l 4 ( Fig. 1c) are different. Indeed, as long as X < 1/2 decreases, the size of a trajectory increases, while its shape changes and looks like a tadpole, with the head centered on L 4 and the tail that extends toward L 3 . In other words, by decreasing X , R oscillates about 1 with an amplitude that increases but remains much smaller than 1, while φ encompasses 60 • with increasing oscillations included in the range ]0 • , 180 • [.
Other families of periodic orbits exist, and several classifications have been realized (see [14,36,37]). Among them, the family f is especially remarkable: it is a one-parameter family of symmetrical periodic orbits whose motion in the synodic reference frame looks like the one of a retrograde satellite of the planet, and that extends from an infinitesimal neighborhood of the planet (i.e., inside its Hill's sphere) to the collision with the Sun (i.e., far beyond the Hill's sphere of the planet). Its computation is similar to the one of the Lyapunov families. However, since it does not originate from a Lagrange fixed point, the initial guess of the method is given by the two following limit cases (see [4]): By varying X > 1, their shape has the same evolution to the one of L 3 . More precisely, the family f seems the symmetrical family of L 3 with respect to the Yaxis, characterized by φ that oscillates about zero and thus a guiding center located on the planet.
The linear stability character of a periodic orbit can be deduced from the monodromy matrix. For ε small enough 2 , the family f is normally elliptic except in two particular orbits that split the neighborhood of the family in three different domains (see [28] for more details). One belongs to the Hill's sphere and corresponds to the "satellized" retrograde satellite orbits. The two others stand for the quasi-satellite orbits, also known as distant retrograde orbits (DRO). Examples of "satellized" retrograde satellite and quasi-satellite orbits, computed during 100 revolutions of Jupiter, are depicted in Fig. 2a-b. Both of the families L l j and L s j are normally elliptic close to the equilibrium. The tadpole-shaped trajectories depicted in Fig. 2c and d, start in the neighborhood of periodic orbits that belongs to L l 4 and L s 4 , respectively. A part of L 3 near the equilibrium is normally hyperbolic and two types of dynamics can be observed in its neighborhood: tadpoleshaped orbits with a large amplitude (Fig. 2e), and horseshoe-shaped orbits that encompass the three fixed points L 3 , L 4 and L 5 (Fig. 2f). The horseshoe-shaped orbit is characterized by R that oscillates about 1 with an amplitude smaller than 1, and φ that features very large oscillations centered on 180 • .
To summarize the situation, we described four types of dynamics-"satellized" retrograde satellite, quasisatellite, tadpole motion and horseshoe motion-that starts in the vicinity of periodic orbits that belong to the six families mentioned above. These dynamics are related by the same features, that is, R that oscillates about 1 and φ that does not circulates but oscillates around a given value. A natural issue is to understand how these dynamics are organized in the phase space of the restricted three-body problem, and especially, if some boundaries can be identified. Nevertheless, the four dimensions of the phase space make difficult the achievement of this goal.
A way to overcome this difficulty is given by a suitable perturbative treatment that focuses on the families of periodic orbits. First of all, let us recall that the particle and the planet are considered in mean-motion resonance, and especially in p:q mean-motion resonance, if they complete, respectively, p and q revolutions around the Sun in the same time. According to the Poincaré classification (see [5]), a periodic orbit of the second or the third "sort" (also translated as "kind") is the continuation, from the limit case ε = 0, of a heliocentric Kepler orbit in mean-motion resonance with the planet. For instance, the families L 3 , L s j as well as the part of the family f that stands outside the Hill's sphere are the continuation of Kepler orbits in 1:1 mean-motion resonance (see the book of Hénon [14] for complete details on the periodic orbit classification). Hence, a perturbative treatment that considers ε as a small parameter, and focuses on a small enough neighborhood of a given mean-motion resonance provides another way to approach some families of periodic orbits and thus to understand the corresponding dynamics. This is the underlying idea associated with the averaged problem that is considered in this work and that we recall in the following section.
2.3 Perturbative treatment of a mean-motion resonance: the averaged problem From now on, we go back to the general case of the restricted three-body problem. If we consider ε as a small parameter, the Hamiltonian function given in the heliocentric reference frame, Eq. (2), can be split in two terms, namely H = H K + H P such that H K corresponds to the unperturbed Kepler motion of the particle, more precisely the motion around a fixed center of mass 1, while H P models the perturbations that depend on ε: the gravitational influence of the planet, the acceleration of the heliocentric frame and a term associated with our choice of Kepler problem. A closed solution of H K describes an ellipse whose shape, orientation and position at a time t are given by the orbital elements (a, e, I, Ω, ω, v(t)). We recall that the position at a time t can also be described by the mean anomaly M(t), a fictitious angle, linear with respect to the time and whose rate of variation-generally known as mean motion-readsṀ(a) = 1/ √ a 3 in the units adopted here. Instead of using the orbital elements, the Poincaré complex variables are adopted because of their regularity for circular and coplanar orbits and since they preserve the symplectic geometry of the problem. In the following, the angles = Ω + ω and λ = M + denote, respectively, the longitude of the periaster and the mean longitude. The symplectic transformation associated with the Poincaré variables readŝ that are, respectively, conjugated to λ,x = −i x and y = −i y. We specify that x and y derive from the angular momentum in the heliocentric reference frame as follows: Moreover, x √ 2/Λ and y √ 8/Λ are equivalent to e exp(i ) and I exp(iΩ) for quasi-circular and quasiplanar orbits.
In the extended phase space, the integrable motion is given by the HamiltonianΞ +Ĥ K withĤ K (Λ) = −1/(2Λ 2 ). Notice that the variable linked to the time is an action variable since the system is periodic. That is the reason why the unperturbed Hamiltonian is integrable. Being the planet and the particle coupled but independent, the solutions describe two ellipses whose respective mean motions are equal toλ(Λ) = 1/Λ 3 anḋ λ = 1. In other words, the solutions of the problem correspond to quasi-periodic orbits with two frequencies.
Since the frequencies are commensurable, the orbits can be periodic. In such a case, the planet and the particle are considered in mean-motion resonance. Studying the restricted three-body problem in this perturbative framework consists in understanding how the pertur-bationĤ P = H P •Υ transforms the unperturbed phase space. More precisely, the analysis can be split in two disjoint situations: -At a suitable distance to mean-motion resonances via a secular model built in order to study the persistence of the quasi-periodic solutions; -Or, on the contrary, in a neighborhood of a meanmotion resonance with the help of special variables and an adapted averaging process.

The resonant variables
For p and q two coprime positive integers, let us consider a neighborhood of the p:q mean-motion resonance. An unperturbed solution is associated with the p:q mean-motion resonance if the semimajor axis of the particle is equal toã = (q/ p) 2/3 . In what follows, a will be defined as the semimajor axis of the "exact" mean-motion resonance. The symplectic transformatioň introduces the resonant angle 3 θ that characterizes the commensurability, and u, its conjugated action, whose modulus measures the distance to the "exact" meanmotion resonance. We recall that θ is not a physical angle and thus it is difficult to represent. Nevertheless, for quasi-circular and quasi-planar orbits, the angular separation φ between the particle and the planet is equivalent to θ + ( p − q)q −1 λ . In the resonant variables given by Eq. (5), the integrable Hamiltonian reads Ξ + H K with H K highlights that θ is constant for u = 0, while it circulates for |u| > 0. More precisely, the angular variables evolve at different rates: λ is a "fast" angle with a frequency 1, θ undergoes "slow" drift in O(u) while ( , Ω) are fixed. Consequently, for a small enough |u|, the timescales of the integrable problem are separated.
In the full problem that reads Ξ + H with all the variables might vary and the motion is very tricky to understand. However, for ε and |u| small enough, the timescales separation between the "fast" and "slow" degrees of freedom still remains. A classical way to exploit this feature is to replace the original problem by another one in which the fast oscillations have been removed. For that purpose, an averaging over the period of revolution of the planet is performed. This process defines the averaged problem.

The averaged problem
The averaged Hamiltonian reads H = H K + H P with choice is arbitrary, and does not affect the dynamics of the solutions but changes their representation. In the framework of a p:q retrograde mean-motion resonance (sometimes denoted as a − p:q mean-motion resonance), the considered resonant angle is neither θ nor proportional to θ (see, e.g. [21,34]). In such a case, that has become an important topic in recent years, the particle orbits around the Sun in the opposite direction to the one of the planet, and the canonical variables usually used in order to describe the motion are not the Poincaré variables introduced by the transformationΥ . However, being the study equivalent to the one of trajectories in (prograde) mean-motion resonance with inclinations |I | > π/2, the following discussions and results remain valid for a retrograde mean-motion resonance. Only the representation of the dynamics will change.
Since H does not depend on λ , Ξ is a first integral that can be dropped. Hence, only three degrees of freedom are required in order to explore the averaged phase space: the resonant variables (θ, u), and (x, x), (ỹ, y), that are, respectively, the Poincaré variables associated with eccentricity and inclination.
There exist at least two classical techniques of computation of the averaged problem. The analytical one is based on the expansion of the Hamiltonian in power series of eccentricity and inclination (see, e.g., [30]). In spite of its efficiency for quasi-circular and quasicoplanar orbits, reaching higher values of eccentricity or inclination requires high-order expansions which generate very heavy expressions. Also worth mentioning the asymmetric expansion developed by Ferraz-Mello and Sato [11] in order to deal with highly eccentric trajectories in mean-motion resonance. The other technique consists on a numerical evaluation of the integral of Eq. (6) and its derivatives. It is a powerful tool since it deals with the Hamiltonian in its exact form which allows to explore the phase space for all values of eccentricity lower than one and all values of inclination.
Following the idea of Poincaré [27], Schubart [31][32][33] developed a numerical averaging procedure for the Hamiltonian in canonical resonant variables. Moons [19] extended the method of Schubart and provided an algorithm that allows to compute the equations of motion of the averaged problem, and thus to construct an integrator for trajectories in p:q mean-motion resonance, for p = q. This algorithm has been adapted by Nesvorný et al. [24] in order to deal with the 1:1 mean-motion resonance. In either case, the numerical averaging can be implemented as follows. Let f an auxiliary function that depends on (θ, u,x, x,ỹ, y, E, E ) where E and E are, respectively, the eccentric anomaly of the particle and the planet. The averaging of f over λ being calculated for fixed (θ, u,x, x,ỹ, y), the Kepler equation implies that: Moreover, since then E can be expressed in terms of (θ, u, x, E) and e . Finally, the averaging of 2π f , that reads is computed by discretizing the variable E as E k = k2π/N with 100 ≤ N ≤ 300 (see [31] for more details).
In the averaged problem, the phase space to explore is 6-dimensional. However, and similarly to the classical approach, the dimension can be reduced in the framework of the circular case (e = 0).

Reduction in the circular case (e = 0)
First of all, we recall that the perturbation H P is analytical outside the collision manifold and thus can be expanded in power series of eccentricity and inclination. In the Poincaré complex variables, the expansion reads where the integers occurring in these summations satisfy the relations In other words, the angular part of the averaged Hamiltonian depends on the linear combination of only two angles: a "modified" resonant angle θ − (q − p)q −1 and the argument of periaster ω = − Ω. Since the averaged Hamiltonian is invariant under the rotations of a third angle, the symplectic geometry imposes the following quantity to be a first integral. We point out that this property can also be derived from the Jacobi constant, Eq. (3).
Indeed, for a given c, such that C (R,Ṙ) = c, the composition of transformations Υ SF •Υ •Υ in resonant variables provides the following expression of the Jacobi constant: Thus, the average of the Jacobi constant over λ introduces the averaged Hamiltonian H , such as c = 2 √ã + ε − 2(H + K ), that is conserved in the averaged problem and implies that K is also a first integral of the averaged problem.
Without revealing details on the symplectic transformation that takes advantage of K , we outline that q( p − q) −1 K , |x| 2 and |y| 2 , respectively, conjugated with θ , − + q(q − p) −1 θ and −Ω + q(q − p) −1 θ, are action-angle variables that can be used. Since these previous variables are singular for the 1:1 mean-motion resonance, the canonical variables that can be adopted are u, K and |y| 2 , respectively, conjugated to θ , − and ω.
Being the degree of freedom associated with K separable to the other two, a reduction is possible. By fixing a value K , seen as a parameter, and eliminating its conjugated cyclic angle, one degree of freedom is removed. Consequently, the averaged phase space can be described by a 1-parameter family of reduced Hamiltonians with two degrees of freedom. We point out that in the circular-planar case, the number of degrees can be reduced to one. Hence, for a fixed K , the "reduced" averaged Hamiltonian is integrable and the description of the phase portrait obtained for various values of K allows to understand the global dynamics of the meanmotion resonance.

Some conclusions about the averaged problem
In Table 1, we summarize the respective features of the averaged problem with respect to the classical approach in the synodic reference frame. First of all, the averaged problem has the advantage to describe the solutions in terms of orbital elements (or variables close to these ones), and thus profits of the symplectic geometry of the problem which allows to reduce by one unit the dimension of the phase space to explore in any case. The algorithms of Moons [19] and Nesvorný et al. [24] are easy to implement to this end. Furthermore, it gives a complete understanding of the resonant dynamics in the circular-planar case. Table 1 Comparison of the features of the restricted three-body problem in the averaged phase space and in the synodic reference frame. C , H and K denote, respectively, the Jacobi constant, the averaged Hamiltonian and the quantity conserved in the aver-aged problem in the circular case. The notations "(N d.o.f)" and "(Non-aut.)" stand, respectively, for "N degrees of freedom" and "non-autonomous" However, the averaged problem possesses also some important drawbacks. First of all, since H has been replaced by H in order to remove the fast oscillations, it does not correspond to the original problem but approximates it with an accuracy that depends on the size of ε. Besides, according to the remark of Schubart [31], it has been shown by Robutel and Pousse [30] and Pousse et al. [28], that the averaged problem fails to describe trajectories that feature close encounters with the planet. In such a case, the "distance" between the averaged Hamiltonian and the original one is important, and the results given by the averaged problem may not be reliable. The clarification of the accuracy as well as the limit of validity is a serious issue. We devote the next section to that purpose.

Notations
Before going further, let us introduce some useful notations associated with the Hamiltonian formalism.
In the following, since it will be necessary to switch from the resonant variables to the heliocentric coordinates, we denote the composition of transformationsΥ •Υ . These two sets of variables preserve the symplectic form, that is, Hence, the Lie derivative of an auxilliary function G(r,ṙ, λ ,Ξ ) along the Hamiltonian flow of a given function F(r,ṙ, λ ,Ξ) reads: t (X 0 ) denotes the Hamiltonian flow at a time t, generated by an auxiliary function h(X) that crosses X 0 at t = 0. We recall that it is connected to the Lie derivative through the exponential of an operator, denoted exp, such as:

The averaging process
According to the perturbation theory, the averaging process coincides with the existence of a symplectic transformation Υ , close to the identity, which maps the original Hamiltonian Ξ + H to Ξ + H + H * , where H * is a remainder that is supposed to be small with respect to H P and thus neglected in the averaged problem. Υ is computed with the time-one map of the Hamiltonian flow generated by some auxiliary function S, that is, In this paper, we choose Based on the previous assumptions, the remainder of the averaging process reads H * can be neglected if and only if it is a perturbation of higher order with respect to H P . However, since H * depends on the derivatives of H P and S that increase as long as the planet and the particle are getting closer, then |H * | and H P can increase simultaneously according to the distance to the singularity and can be at least of the same order. In such a case, the hierarchy between the perturbative terms is not ensured and the approximation provided by the averaged Hamiltonian H might not reflect properly the dynamics of the restricted three-body problem. In other words, in the neighborhood of the collision manifold can exist an "exclusion zone" inside which the solutions of the restricted three-body problem fall outside the scope of the averaged Hamiltonian.
The following section is devoted to the characterization of this exclusion zone through a quantitative treatment of the averaging process.

Quantitative treatment of the averaging process
We first introduce a domain and a norm on the extended phase space that will allow us to compute quantitative estimates.
In this development, we are not interested in the situations of close encounters with the Sun, that occur for a very high eccentricity. They are avoided by consideringΔ and σ as arbitrarily fixed small numbers independent of ε, ρ and Δ.
The estimates will be computed through the supremum norm on D κ , denoted · κ such that i≤n is a n-dimensional vector field that depends on the resonant variables (θ, u,x, x,ỹ, y, λ ,  Ξ). Since we do not attempt to obtain estimates with particularly sharp constants, all constants have been suppressed and replaced by the Pöschel's notation, that is, x ≤ • y, x • ≤ y, and x = • y to indicate, respectively, that x < Cy, C x < y, and x = Cy with some constant C ≥ 1 independent of ε, ρ and Δ.
In this setting, the size of the functions involved in the averaging process can be estimated. Hence, we state the following: Consequently, H K , H P , H P and Υ are bounded together with their partial derivatives with respect to θ , u,x, x,ỹ and y. More precisely, for each order k lower or equal to an arbitrarily fixed n ≥ 1 and (W i ) i≤k ∈ {θ, u,x, x,ỹ, y}, the following thresholds are satisfied on D 2 : The previous lemma allows to state an averaging theorem where quantitative estimates on the averaging process are computed.

Theorem 1 Assuming ε, ρ and Δ small enough such that
there exists a symplectic transformation close to the identity, denoted as Υ : , u . ,x . , x . ,ỹ . , y . , λ , Ξ . ) → (θ, u,x, x,ỹ, y, λ , Ξ), for W ∈ {θ, u,x, x,ỹ, y}, such that, in the "averaged" resonant variables (θ . , u . ,x . , x . ,ỹ . , y . , λ , Ξ . ), the Hamiltonian reads: where H * is the remainder of the averaging process. Furthermore, on the domain D 3/2 , H * together with its partial derivatives with respect to θ , u,x, x,ỹ and y are bounded and satisfy the following thresholds: Lemma 1 and Theorem 1 provide the estimates that allow to compare how H P and |H * | increase as long as the planet and the particle are getting closer. In order to understand the role of the distance to the planet in the size of the remainder, we first relate the upper bound of the distance to the resonance to ε and Δ by choosing which satisfies Eq. (9) and such that the two terms in η depend on the same quantity: the ratio ε Δ 3 . Hence, we have which imposes the lower bound Δ ≥ • ε 1/3 in order to get decreasing perturbations in the "averaged" resonant variables. More precisely, we recover the size of the Hill's sphere of the planet. We recall that Δ denotes the minimal mutual distance, that is, the minimal distance between the particle and the planet, which is allowed in the considered domain D κ . As a consequence, if we relate Δ to ε as follows: 1/3 the Hill's radius of the planet, then Theorem 1 ensures that the domain D 4/3 stands outside the exclusion zone of the averaged problem for ε small enough and N ε > 1, that is, for a minimal mutual distance larger than the Hill's radius of the planet. In spite of this feature, Theorem 1 does not establish that a given solution of the averaged problem that starts inside D 3/2 , does not escape and does not cross the exclusion zone of the averaged phase space. To this end, a careful analyze of the behavior of the averaged solutions has to be led for each type of dynamics in mean-motion resonance. Nevertheless, assuming that the solution remains inside D 3/2 until a certain amount of time, a theorem of stability over finite times can be proven in the restricted three-body problem.
Before stating the theorem, let us denote a solution governed by the averaged Hamiltonian Ξ + H , that starts in X 0 ∈ D 1 and remains inside this domain up to a given time T 1 > 0, as and λ (t) = t.
For |t| ≤ T 1 and N ε > 1, X . (t) approximates the solution of the restricted three-body problem that starts in X 0 . In the resonant variables, the solution governed by Ξ + H can be written as (W(t), λ (t), Ξ (t)) with being (δ i (t)) i≤6 the functions that denote the error on each coordinate in the approximate solution given by the averaged problem.

Theorem 2
With the previous notations, the errors in the approximate solution satisfy the following upper bounds: for |t| ≤ min (T , T 1 ) and T = 2π N 3 ε .
In the following reasonings, it is assumed that T ≤ T 1 . First of all, in the limit case given by N ε = ε −1/3 , that is, Δ = 1 and ρ = • √ ε, Theorem 2 ensures that, up to a finite time in O(1/ √ ε), the approximate solution given by the averaged problem remains inside a neighborhood of order O( √ ε) of the solution obtained in the synodic reference frame. Hence, we obtain the results for which the particle is considered distant enough from the planet. The other limit case, that is, for N ε 1, the particle can approach the edge of the Hill's sphere while the distance to the resonance can reach the order O(ε 1/3 ). Even though the gravitational influence of the planet is not dominant, it can be strong enough with respect to the one of the Sun. Hence, we can only ensure that the accuracy of the approximate solution will not exceed a quantity of order O(ε 1/3 ) for one or few periods of revolution of the planet. In such a case, the solution of the averaged problem might not be reliable in order to approach the one obtained in the synodic reference frame.
By increasing N ε > 1, the minimal mutual distance moves away from the Hill's sphere, the gravitational effect of the planet becomes weaker with respect to the one of the Sun, and the upper bound on the distance to the mean-motion resonance decreases as ρ = • ε 1/3 √ N ε . Thus, the approximate solution becomes more accurate with an upper bound on the error that decreases as 1/ √ N ε for the resonant angle θ and 1/N 2 ε for the other variables, and a time of stability that increases as N 3 ε . Consequently, multiplying by a factor 5 the minimal mutual distance divides the errors by a factor √ 5 2 and 25, respectively, on θ and the other variables, and multiplies the time of stability by a factor 5 3/2 10.
In a more practical way, for a given number n > 0 of revolutions of the planet, with n ≤ O(1/ √ ε), this result provides a domain of initial conditions such that the approximate solution given by the averaged problem is reliable in order to approximate the one obtained in the synodic reference frame. The errors in the approximate solution being, respectively, ε 1/3 n −1/3 in θ and ε 1/3 n −4/3 in the other variables, the smaller the planet mass ratio is, the greater the accuracy would be.

Discussion
The proofs of Lemma 1, Theorems 1 and 2 are given in the Appendix A.
The key ingredient of the proof of Lemma 1 is our definition of the collisionless domains D κ given in terms of heliocentric coordinates instead of resonant variables in order to exclude a neighborhood of the collision manifolds. Since the Poincaré complex variables prevent singularities associated with the eccentricity or inclination equals to zero, the transformation Υ is analytic and can be bounded independently to ε, ρ and Δ. Hence, the estimates on H K , H P and H P are deduced directly from H K and H P .
The proof of Theorem 1 has two parts. We first characterize the conditions that allow to ensure that the transformation of averaging Υ is close to identity and maps the domain D 3/2 in the domain D 2 inside which the estimates are computed in Lemma 1. In the second part, we estimate the remainder H * , Eq. (8), and its associated vector field by using the Taylor expansions at zero and first order combined with the estimates of Lemma 1.
Finally, Theorem 2 is a direct application of the classical strategy to prove stability over finite times (see [2]). For that purpose, we compare the vector field of the approximation given by the averaged problem with the one of a solution of the original problem. Assuming that the two solutions remain in a given neighborhood up to a time T , we can choose T such that the order on the errors in the approximation is of the same order as the one of the transformation in "averaged" resonant variables.
The validity of the averaged problem was discussed by Robutel and Pousse [30], Robutel et al. [29] and Pousse et al. [28] in the framework of the 1:1 meanmotion resonance. Moreover, it is a key ingredient of the proof given by Niederman et al. [25] on the existence 4 of the horseshoe-shaped trajectories followed by the two Saturn's moons, Janus and Epimetheus. Since the semimajor axes of the two small bodies are almost the same, the issue generated by periodical close encounters is manifest. For instance, in [30] and [28], the authors highlight, through a frequency analysis, that, when an initial condition tends to the singularity of mutual collision, the approximation given by the averaged problem has fundamental frequencies that increase and tend to infinity. This phenomenon is inconsistent with the hypothesis of timescales separation required by the averaged problem. For that reason, an arbitrary criterion 5 on the frequencies was introduced in [28] in order to localize the exclusion zone in the averaged phase space.
The work of Robutel et al. [29] provides a rigorous treatment of the averaging process for resonant dynamics in the planar planetary three-body problem. In the framework of quasi-circular co-orbital trajectories, it gives quantitative estimates on the remainder generated by the averaging process as well as its vector fields in order to state a theorem of stability over finite times. The proof uses a complex domain of holomorphy inside which estimates are computed through the Cauchy inequality. The results of the present paper are based on the same idea but applied in the general case (eccentric and spatial trajectories) of a generic p:q mean-motion resonance of the restricted three-body problem. Even though the technique of complexifying used in [29] is very efficient in the case of quasicircular and quasi-planar orbits, the definition of the minimal mutual distance Δ in terms of resonant variables is really difficult in the general case. That is why, we have chosen the direct computation of estimates by taking advantage of the form of H K , as well as the one of the H P , that only depends on r and λ .

The co-orbital motion in the circular-planar case
In this section, we focus on the co-orbital motion (1:1 mean-motion resonance) in the circular-planar case. More precisely, in the framework of the averaged problem, we intend to approach the six families of periodic orbits described in Sect. 2.2.2 (the short-periodic L s j and long-periodic L l j for j = 4, 5, the Lyapunov family L 3 , and the family f ) as well as the dynamics observed in their neighborhood (the tadpole motion, the horseshoe motion, the quasi-satellite motion, and the "satellized" retrograde satellite motion), and identify the limit of validity of the corresponding solutions by applying Theorems 1 and 2 developed in Sect. 3.
First of all, we introduce the resonant variables and apply the properties stated in Sect. 2.3 to the case of the 1:1 mean-motion resonance. The resonant degree of freedom is described by the angle θ = λ − λ and the action u = √ a − 1 that measures the distance to the 5 A solution of the averaged problem was considered outside the exclusion zone, if the modulus of their fundamental frequencies are lower thanλ /4 whereλ = 1 denotes the frequency of averaging. exact mean-motion resonance given by the semimajor axisã = 1. In the circular-planar case, the secular variations of the orbits are provided by (x, x) for which K = |x| 2 is a first integral of the averaged problem. For a fixed K , seen as a parameter, the reduced averaged Hamiltonian, denoted as follows is integrable with one degree of freedom, and thus allows to understand the co-orbital dynamics through a phase portrait. Instead of using K , we introduce the parameter e 0 such as Hence, e 0 defines the eccentricity of a trajectory that crosses the orbit of the planet, that is, at the exact meanmotion resonance u = 0. Furthermore, since e 0 approximates the eccentricity of the trajectories that belong to a given phase portrait. Notice that e 0 is also connected to the Jacobi constant through Eq. (7). More precisely, Theorem 1 and Eq. (11) imply that inside the domain of validity of the averaged problem, denoted D 1 (Δ) in the previous section, the Jacobi constant C (R,Ṙ) = c reads: Consequently, Eq. (10) and the estimates of Lemma 1 ensure the following relation between the two quantities: Before going further, we will see in the next section how a trajectory of a given phase portrait is related to the solutions of the averaged problem. Besides, in order to bridge a gap between the classical and the perturbative approaches, we will detail how a solution of the averaged problem describes the motion of a particle in the synodic reference frame. These correspondences were described in [28], but the relationship with the synodic reference frame needed to understand the dynamics were lacking.

Reading a phase portrait
For a given value of e 0 ≥ 0 and a given initial condition (θ i , u i ), a trajectory that belongs to the corresponding phase portrait is generally 6 a periodic solution but can also be an equilibrium of H K (e 0 ) . If we denote its frequency ν, a periodic trajectory can be written as where the functions F j are 2π -periodic such that F j (0) = 0. The dynamics of the angle = arg(x) is required in order to relate this trajectory to the solutions of the averaged problem. Sincė is 2π/ν-periodic, there exists F 3 , a 2π -periodic function with mean zero and F 3 (0) = 0, such that for all is the secular precession frequency of . In other words, a solution of the averaged problem that starts can generally be written as and a periodic trajectory of a given phase portrait generally corresponds to a set of quasi-periodic solutions parametrized by i ∈ T, whose fundamental frequencies are given by ν and g. The same reasoning applies for an equilibrium of H K (e 0 ) : it corresponds to a set of periodic solutions of the averaged problem, parametrized by i ∈ T, and that can be written as Nevertheless, being ignorable when the osculating ellipses are circles (e 0 = 0), the solutions have the same features in the averaged problem as in the phase portrait of H K =0 .
Theorem 2 ensures that a given solution of the averaged problem that lies outside the Hill's sphere of the planet, approximates for a finite time the motion of a particle that starts at the same initial condition. Hence, in the Poincaré complex variables, the motion of a particle that crosses (λ i , Λ i ,x i , x i ) with λ i = θ i and Λ i = 1 + u i at t = 0, can be approximated for a finite time by Eq. (12) or Eq. (13) In terms of orbital elements, the approximation of the variations reads Hence, the semimajor axis and the eccentricity experience a slow oscillation of frequency ν and of amplitude of the order O(u), respectively, about 1 and e 0 . The variations of the longitude of the periaster correspond to the composition of a secular drift of frequency g with an oscillation of frequency ν. Finally, the motion of the mean longitude is given by a fast drift of frequency 1 − g composed with a slow oscillation of frequency ν.
We recall that the orbital elements are related to the polar coordinates (φ, R) = (arg(R), R ) of the synodic reference frame as follows: M). (14) The functions G j , that satisfy G j (e, M) = O(e), derive from the Kepler equation M = E − e sin E and the difference between the true and eccentric anomaly tan(v/2) = 1+e 1−e tan(E/2). In the synodic reference frame, the approximate motion of the particle can be written as As a consequence, a periodic trajectory of a given phase portrait generally provides a quasi-periodic approximation of the motion whose fundamental frequencies are 1 − g and ν. More precisely, the motion is characterized by R that oscillates about 1 with an amplitude of the order O(e 0 ) + O(u(t)), while φ is the sum of the periodic oscillation generated by θ(t) with the quasiperiodic oscillations generated by G 1 whose amplitude is of the order O(e 0 ). The same reasoning applies for an equilibrium of H K (e 0 ) : it provides a set of periodic trajectories of frequency 1 − g whose motion follows: Hence, the motion is characterized by an oscillation of frequency 1 − g around a guiding center located in (φ, R) = (θ i , 1 + O(u i )) and whose amplitude is of the order O(e 0 ). Finally, for e 0 = 0, the motion of the particle is approximated by and a trajectory of the phase portrait provides an equilibrium in the synodic reference frame or a periodic trajectory of frequency ν characterized by R that oscillates around 1 with an amplitude of the order O(u), that is, much smaller than 1.  [24]. They are equivalent to the one showed in [22,24] and extensively described in [28]. We will limit ourselves to present what will be useful to understand the trajectories described in Sect. 2.2.2 as well as to apply the Theorems stated in Sect. 3. Notice that the phase portraits are invariant by the symmetry with respect to the u-axis (θ, u) → (2π − θ, u).
In Fig. 3, e 0 is equal to zero and particle's motion is quasi-circular. Let us mention that in this case the reduced averaged Hamiltonian reads 2 , which allows to recover some classical analytical properties (see, e.g., [12,23]) that will be recalled in the following. Being the computations similar to the ones given in [30], the details are not given.
First of all, the black dot at θ = u = 0 embodies the collision with the planet, where H 0 is not defined. The phase portrait possesses five equilibria that correspond to the Lagrange fixed points L j . The two elliptic equilibria located in (θ j , u j ) = (−1) j × 60 • , 0 , stand for L 4 and L 5 while the hyperbolic ones in The phase portrait is characterized by six regions. The beige domain centered at the singularity and bounded by the separatrices originating from L 1 and L 2 seems to be the prograde satellite-like motion. However, since the domain belongs to the Hill's sphere of the planet, it lies inside the exclusion zone and its trajectories fall outside the scope of the averaged problem. Hence, the corresponding dynamics will not be analyzed in this work. The upper and lower grey domains, that lay above the separatrices of L 1 and L 2 , illustrate the non-resonant motion for which θ circulates (clockwise in the upper region and anti-clockwise in the lower one). Considering that the width along the u-axis of the region located inside the separatrices is of the order O(ε 1/3 ), Theorem 1 implies that the nonresonant regions escape from the domain of validity of the averaged problem.
The three remaining regions are the ones of the coorbital dynamics, for which θ oscillates about a given Fig. 3 Phase portrait of a particle in quasi-circular motion (e 0 = 0) for a Sun-Jupiter-like system (ε = 1/1000). The black dot denotes the singularity associated with the collision with Jupiter, while the brown circle, orange circle, red circle and the two blue diamonds correspond, respectively, to L 1 , L 2 , L 3 and L j for j = 4, 5. The separatrices that originate from L 1 , L 2 and L 3 , represented, respectively, by brown, orange and red thick curves, divide the phase portrait in six regions. Although they are very close to each other, the separatrices of L 1 do not coincide with those that originate from L 2 . The beige domain, centered on the singularity, embodies the Hill's sphere of the planet, inside which the averaged Hamiltonian does not reflect properly the dynamics of the restricted three-body problem. The upper and lower grey domains lay outside the co-orbital resonance (θ circulates). The blue and red trajectories are level curves associated with tadpole-shaped and horseshoe-shaped periodic orbits, respectively. More precisely, the blue level curves belong to the families L l 4 and L l 5 . Finally, the green curves exhibit the elements associated with a minimal mutual distance Δ = N R H , while the grey lines depict {|u| = ε 1/3 / √ N }. For a given N , they bound the domain D 1 (Δ) inside which Theorem 2 is applied . (Color figure online) value. These regions are divided by the separatrix that originates from L 3 . The solutions that librate around L 4 and L 5 are tadpole-shaped characterized by |θ | that oscillates about 60 • such that Θ 0 < |θ | < 180 • with and u that oscillates around zero with an amplitude that can reach O( √ ε). According to Sect. 4.1, these trajectories are periodic and possess the same features as the ones of the long-periodic families L l 4 and L l 5 . Outside the separatrix, the trajectories encompass L 3 , L 4 and L 5 , i.e., they are characterized by θ and u that oscillate, respectively, about 180 • and 0, with large amplitudes. More precisely, θ and u undergo large variations such These solutions approximate periodic horseshoe-shaped trajectories (see, e.g., the study of Barrabés and Ollé [3] that focuses on these solutions).
For e 0 > 0, the location of the singularities evolves: the origin becomes a regular point surrounded by a set of singular points that describes a curve. For small e 0 (e.g., Fig. 4a), a new domain of co-orbital motion appears inside the collision curve. It is centered on an elliptic equilibrium point located close to the origin. According to Sect. 2.2.2 and Sect. 4.1, the elliptic equilibrium point approximates a periodic orbit of the family f . Hence, the periodic trajectories that librate around, provide quasi-periodic approximations of quasi-satellite orbits. Outside the collision curve, the topology does not change with respect to the one depicted in Fig. 3 outside the Hill's sphere: two elliptic equilibria close to L 4 's and L 5 's locations and a separatrix emerging from a hyperbolic equilibrium close to L 3 that divides the regions of tadpole and horseshoe motions. According to Sect. 4.1, the equilibria of the phase portraits approximate periodic trajectories whose guiding center are located, respectively, close to  Fig. 4. For a, b, c and d, e 0 is equal to 0.15, 0.5, 0.75 and 0.925, respectively. The black curves represent collision with the planet. The blue, sky blue and red trajectories are level curves associated with tadpole, quasi-satellite and horseshoe motion, respectively. The red circles, sky blue and blue diamonds are equilibria. More precisely, they corre-spond, respectively, to orbits of the families of periodic orbits L 3 , f and L s j for j = 4, 5. From each hyperbolic equilibrium emerges a separatrix (red thick curve) that divides tadpole motion and horseshoe motion. For d, the red diamond denotes an orbit that belong to the stable part of L 3 . Finally, the green curves exhibit the elements associated with a minimal mutual distance Δ = N R H , and bound the domain D 1 (Δ) inside which Theorem 2 is applied. (Color figure online) For higher values of e 0 (e.g., Fig. 4b-c), the size of the quasi-satellite domain increases while one of the tadpole domains shrinks when the two elliptic equilibria are getting closer to the hyperbolic one.
The evolution of the two elliptic equilibria illustrates the shift of the guiding center of L s 4 and L s 5 toward Finally, for very high e 0 (e.g. Fig. 4d), the tadpole domains vanished and remains a domain characterized by trajectories that librate around an elliptic equilibrium that belongs to L 3 . This bifurcation of L 3 occurs for e 0 0.917, when the short-periodic families L s j merge with L 3 . For a given minimal mutual distance Δ = N R H with 1 ≤ N ≤ ε −1/3 and R H = ε 3 1/3 that denotes the Hill's radius of the planet, we recall that the domain D 1 (Δ) introduces in Sect. 3.3 is defined by the resonant variables for which the distance between the planet and the particle is larger than N R H , and |u| which is bounded by a quantity ρ = • ε 1/3 √ N . The grey lines and green curves that lay over the phase portraits approximate the boundaries of D 1 (Δ) for several value of N = 1, 3, 5, 10. The grey lines correspond to {|u| = ε 1/3 √ N } while the green curves exhibit the elements (θ, u, e(e 0 , u)) for which the minimal mutual distance is equal to Δ = N R H . They are computed by resolving the following equation: with the help of Eq. (14).
In Figs. 3 and 4, the exclusion zone of the averaged problem is depicted by the small areas centered on the collision curves and bounded by the continuous green curves associated with the minimal mutual distance Δ = R H . They show that, contrarily to the tadpole motion, some solutions in quasi-satellite and horsehoe motion intersects the Hill's sphere and fall outside the scope of the averaged Hamiltonian. More generally, the phase portraits also show that, for a given minimal mutual distance Δ = N R H , a co-orbital solution which starts inside the area {|u| ≤ ε 1/3 √ N } can cross the sphere defined by Δ = N R H and thus escape from the domain D 1 (Δ) inside which Theorems 1 and 2 are applied. However, a co-orbital solution which starts in the neighborhood of the section {u = 0} at a given minimal mutual distance Δ, does not experience closest encounters with the Jupiter and remains inside D 1 (Δ). As a consequence, the section {u = 0} provides a convenient way to discuss about the validity and the time of stability of the solutions of the averaged problem without caring about times of escape from D 1 (Δ). This study is realized in the following section.

A "map" of the co-orbital motion in the circular-planar case
First of all, we summarize the situation. Six families of periodic orbits ( f , L 3 , L s j and L l j ) and three types of trajectories (tadpole, horseshoe and quasi-satellite motion) described in Sect. 2.2.2 have been recovered in the averaged problem close to the exact mean-motion resonance u = 0. Notice that the domain of "satellized" retrograde satellite mentioned in Sect. 2.2.2 is missing since it is located in the neighborhood of the family f that belongs to the Hill's sphere (see [28] for more details).
Each domain of co-orbital motion extends quasisymmetrically with respect to the exact mean-motion resonance u = 0 and is neatly defined by the collision curves or the separatrices that originate from the hyperbolic equilibria associated with L 3 . In this section, we construct a "map" of the co-orbital motion in the circular-case, that is, a representation of the section {u = 0} which can be used in order to discuss about the stability of the solutions as well as to compute coorbital trajectories in the synodic reference frame. Two parameters are required to identify a solution of the averaged problem that belongs to the section {u = 0}: the resonant angle θ i and the eccentricity of the orbit e i , (we recall that e i = e 0 when u = 0). Hence, we compute the evolution of the cross sections of the boundaries of each domain (separatrix and collision curves) by varying the eccentricity e i . To that end, we consider the following reduction of the reduced averaged Hamiltonian: which is derived from the Taylor expansions of H P and H K , respectively, at zero and second order in u = 0. This reduction is reliable in the vicinity of the section {u = 0} and has the advantage to be independent of ε, under the following rescaling: Hence, a "map" of the section {u = 0} computed through H K (e 0 ) 0 , is invariant under the variation of ε. Figure 5 displays the map of the co-orbital motion in the circular-planar case. The black thick curves illustrate the collision with the planet. They can be approximated by |θ i | = 2e i × 180 • /π up to high eccentricities. They bound the sky blue region associated with the quasi-satellite motion. The red thick curves depict the crossings of the separatrices that originate  Fig. 6 displays the map of the co-orbital motion with the elements (θ i , e i ) for which the minimal mutual distance Δ is equal to N R H , for a Sun-Jupiter-like system (ε = 1/1000). The two additional grey lines represent the elements (θ i , e i ) for which the particle crosses the orbit of Mars and Saturn. Hence, they suggest the maximal value of eccentricity for which the solutions of the restricted three-body problem are reliable in order to describe the real motion in the Solar System. The right panel of Fig. 6 is an enlargement of the map on the region that surround the collision curve.
According to Theorems 1 and 2, for a given number 1 < N ≤ ε −1/3 , and a given type of co-orbital dynamics, we can define a set of elements (θ i , e i ) which satisfy a mutual distance greater than Δ = N R H , and for which the time of stability of the solutions of the averaged problem is at least √ N 3 revolutions of Jupiter.
In the synodic reference frame, the couple (θ i , e i ) provides a set of initial conditions, parametrized by i ∈ T and that can be written as Hence, Theorem 2 ensures that an initial condition, given by Eq. (16), provides a co-orbital trajectory of the same type, at least for a finite time. In other words, transitions to another co-orbital motion or escapes from the 1:1 mean-motion resonance cannot occur at least during a time T = 2π √ N 3 . For instance, for N equal to 3, 5 and 10, it ensures a time of stability corresponding approximately to 5, 10 and 30 revolutions of Jupiter, that is, more than 50, 125 and 350 years.
As mentioned at the end of the previous section, the quasi-satellite and the horseshoe domains intersect the exclusion zone of the averaged problem. More precisely, for high values of e i , the quasi-satellite motion dominates the map and the size of the intersection between the quasi-satellite domain and the exclusion zone is small relatively to the whole domain. By decreasing e i , since the quasi-satellite domain shrinks with the collision curve, the relative size of the inter-section increases until a critical value for which the exclusion zone contains all the quasi-satellite orbits. In the case of a Sun-Jupiter-like system, this critical value occurs for e i 0.07.
Notice that Pousse et al. [28] suggested, through a frequency analysis of the family f , a critical value of e i 0.18 for the quasi-satellite orbits. This value was given by an arbitrary criterion for which a solution of the averaged problem is considered outside the exclusion zone, if the modulus of their fundamental frequencies ν and g are lower thanλ /4, whereλ = 1 denotes the frequency of averaging. Theorem 1 provides a lower value of eccentricity and thus a larger domain of validity of the averaged Hamiltonian for the quasi-satellite motion. However, based on the results in [28], the quasi-satellite region that surround the Hill's sphere is probably overlaped by secondary resonances, that is, the resonant structures generated by commensurabilities between frequencies ν, g andλ = 1, and especially between ν and 1 − g due to the D'Alembert rules (see Sect. 2.3.2). In particular, it has been shown in this area (see, e.g., [28]) that the neighborhood of the family f is divided in three disjoint regions by two critical orbits of the family f associated with the commensurabilty 3ν = 1 − g. As a consequence, a global study of the frequencies in the averaged problem will probably reveal the resonant structures that destabilizes the quasi-satellite region that surround the Hill's sphere, and may also highlight some islands of quasi-satellite solutions for which the time of stability is larger than the one given by Theorem 2.
The horseshoe motion exists outside the exclusion zone for low and very high eccentricities. More precisely, by increasing e i , the size of the intersection between the horseshoe domain and the exclusion zone increases until a critical value e i 0.4 for which all the horseshoe-shaped trajectories cross the section {u = 0} inside the Hill's sphere. Furthermore, there exists another critical value e i 0.6 for which a part of the horseshoe domain goes outside the exclusion zone. Then, for increasing e i , the relative size of the intersection between the horseshoe domain and the exclusion zone decreases.
Most of the solutions in horseshoe motion experiences closed encounters with Jupiter (less than 5 Hill's radii), and thus have a relatively small time of stability according to Theorem 2. Similarly to the quasi-satellite region located at the edge of the Hill's sphere, the horseshoe region is probably overlaped by secondary res-onances which destabilize the domain. A global frequency analysis of the region may reveal these resonant structures.

Conclusions
In this paper, we showed that the averaged problem provides another approach in order to study some families of periodic orbits of the restricted three-body problem. More precisely, we proved that it is a valid approximation in a particular area of the phase space that focuses on mean-motion resonances. Through a rigorous treatment, we characterized the domain of validity of the averaged problem (Theorem 1) and proved that it is a reliable approximation as long as the considered trajectories lay outside the Hill's sphere of the planet. A new result of stability over finite times has also been proved (Theorem 2). As a consequence, we provided a rigorous justification of the relevance of the averaged problem to study some specific solutions of the restricted three-body problem.
Our theoretical results allowed us to understand the co-orbital motion (1:1 mean-motion resonance), that is, the quasi-satellite, the tadpole and the horseshoe orbits, that comprise the family f , the short-periodic and longperiodic families that originate from L 4 and L 5 , and the Lyapunov family associated with L 3 . In particular, in the framework of the circular-planar case, we propose a method, illustrated by a "map" of the co-orbital motion, that takes advantage of the averaged problem in order to compute co-orbital trajectories in the synodic reference frame. The results are presented in the case of a Sun-Jupiter-like system, but the "map" of the co-orbital motion, plotted in the Fig. 5, is independent of the small parameter ε, that is, the mass ratio of the planet over the total masses of the system. Hence, only the elements for which the minimal mutual distance Δ is equal to N R H , must be computed in order to apply the method to a different Sun-planet system (or planet-moon system). For example, Fig. 7 displays the map of the co-orbital motion for a Sun-Earth like system (ε = 1/333333). However, we recall that, since the accuracy of the averaged problem depends on ε, the larger ε is, the less the map of Fig. 5 is reliable.
A practical application of our theoretical results interests the design of space missions. Let us consider a spacecraft affected by the gravitational forces of a Sun-Earth like system. The selected orbits are usually remarkable solutions in the synodic reference frame that is, the Lagrange fixed points and periodic or quasiperiodic trajectories, as well as their associated hyperbolic manifolds (when existing). However, except for the dynamics associated with L 1 and L 2 , most of these solutions are located at a remote distance from the Earth, and thus the cost in terms of energy or time to reach them is usually not affordable. For instance, Fig. 7 shows that L 3 , L 4 , L 5 , the short-periodic and long-periodic families and the Lyapunov family L 3 lay at a distance larger than 40 Hill's radii. Only the family f provides periodic orbits available at a lower distance, that is why it becomes an important topic for mission design, mainly in the aim to orbit small moons (e.g., the two Mars's moons, Phobos and Deimos).
Our idea is the following: since the duration of a mission is limited, it is not necessary to target an equilibrium or a periodic solution in the synodic reference frame. With the help of our method, which characterizes the elliptic and hyperbolic dynamics of the coorbital motion through the averaged problem as well as defines a time of stability of these solutions, it is easy to select an initial condition on the map of Fig. 7, that is located close enough to the Earth, and that satisfies a given co-orbital dynamics for the whole duration of the mission. For instance, for a 30-years mission, Theorem 2 ensures that a co-orbital solution, that is located outside 10 Hill's radii of the Earth, will be stable at least during the time of the mission. As a consequence, the horseshoe, tadpole and quasi-satellite solutions become possible target trajectories.
Notice that co-orbital solutions that experience closest approaches with the Earth may also be stable for the considered time of the mission. For instance, in the case of planar and quasi-circular horseshoe orbits, Cors and Hall [7] provided a practical stability estimate (the order of magnitude of the time for which the trajectory remains qualitatively the same) which is much greater than the one given by Theorem 2. Hence, a complementary analyze could be a global numerical investigation of the time of stability of the co-orbital orbits that experience close approaches. A detailed study will be addressed in a forthcoming work.
Finally, we point out that the application of the averaged problem as presented here is not restricted to the co-orbital motion in the circular-planar case: it can also be applied to inclined co-orbital trajectories, or, more generally, to solutions associated with other Fig. 7 Map of the co-orbital motion for a Sun-Earth like system mean-motion resonances. To this aim, a careful study of the averaged phase space must be realized.
Funding This work was supported by the project entitled "coorbital motion and three-body regimes in the solar system," funded by Fondazione Cariplo through the program: "Promozione dell'attrattività e competitività dei ricercatori su strumenti dell'European Research Council-Sottomisura rafforzamento" Data availability Data sharing not applicable to this article as no datasets were generated or analyzed during the current study.

Conflict of interest
The authors declare that they have no conflict of interest that are relevant to the content of this article.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/ by/4.0/.

A.1 Remainders
We recall that the Hamiltonian flow at a time t, generated by an auxiliary function h(X), satisfies the following property: where g is an auxiliary function. Thus, we have the following Taylor expansions: A.2 Proof of Lemma 1 For given ρ • ≤ 1, Δ • ≤ 1, ε • ≤ 1 and κ > 0, the domain D κ has been designed in order to exclude the collision manifold and make the perturbation H P analytic. Hence, H P is bounded on D κ as well as its partial derivatives with respect to θ , u,x, x,ỹ or y for each order k lower or equal to an arbitrarily fixed n ≥ 1. Their estimates are deduced from the following reasonings.
First of all, we recall the thresholds provided by the definition of D κ : Δ/κ< r • Υ κ ≤ • 1, Δ/κ≤ r • Υ − r κ ≤ • 1, r κ = • 1, whereΔ = O(1) since it is an arbitrarily fixed quantity that does not depend on ε, Δ and ρ. For all order k ≤ n, the perturbation H P in heliocentric cartesian coordinates, Eq. (4), yields the following estimates on D κ : Since the transformationΥ , that introduces the Poincaré complex variables, Eq. (5), is regular when eccentricity and inclination tend to zero, it does not have singularities. More precisely, it is an analytic transformation on D κ and, for all order k ≤ n, its derivatives can be bounded by a constant that does not depend on ε, ρ and Δ.Υ , which introduces the resonant variables, Eq. (5), is an affine transformation that fulfills the same properties. As a consequence, Υ is analytic and, for each order k ≤ n and (W i ) i≤k ∈ {θ, u,x, x,ỹ, y}, the following threshold is satisfied: Finally, for a given analytic function F(r, λ ) and W ∈ {θ, u,x, x,ỹ, y}, we recall the chain rule: The bounds on the partial derivatives of H P with respect to θ , u,x, x,ỹ and y, and for all order k ≤ n are deduced from the combination of the chain rule, Eq. (22), with the thresholds given by Eqs. (20) and (21).
Since H P κ ≤ H P κ , the results on the averaged perturbation is a direct consequence of the previous development.
On the domain D 2 , H K is analytic, is different from zero, and thus satisfies H K 2 = • 1 while, for each order k lower or equal to an arbitrarily fixed n ≥ 1, its derivatives are bounded. More precisely, in order to satisfy the following property: Under these conditions, the remainder of the averaging process reads: Most of the estimates required for the proof are computed with the Taylor expansion at zero order, Eq. (18). More precisely, if we assume that there exists κ < 2 such that Υ (D κ ) ⊂ D 2 , then for a given function g that depends on the variables (θ, u,x, x,ỹ, y, λ ), the following thresholds are ensured: For each order k lower or equal to an arbitrarily fixed n ≥ 1 and (W i ) i≤k ∈ {θ, u,x, x,ỹ, y}, we point out that Eq. (23) and Lemma 1 provide the following thresholds on the partial derivatives of S: Hence, in D κ and for W ∈ {θ, u,x, x,ỹ, y}, we have the following: As a consequence, for ε • ≤ ρΔ 2 and ε • ≤ Δ 3 with ε, ρ and Δ small enough, we can choose κ = 3/2 such that the symplectic transformation of averaging satisfies Υ (D κ ) ⊂ D 2 and is close to identity such that for W ∈ {θ, u,x, x,ỹ, y}. It remains to prove the estimates on the remainder of the averaging process. The Taylor expansions at zero and first order, Eqs. (18)(19), combined with the condition satisfied by S, Eq. (24), ensure that the remainder (1 − s)L S H P • Φ S s ds.
As a consequence, the thresholds given by Eq. (26) and Lemma 1 provide the following upper bound on the remainder: The upper bound on the derivative of H * with respect to θ , u,x, x,ỹ and y is deduced in the same way. For W ∈ {θ, u,x, x,ỹ, y}, Lemma 1 and Eqs. (26)(27) provide the following estimates: x, x,ỹ, y}.

A.4 Proof of Theorem 2
The proof follows the classical strategy described in [2]. The aim of the first part of the proof is to bound the Hamiltonian vector field linked to H * over the domain D 3/2 . The second part comes from the size of the transformation of averaging Υ and the choice of a time T > 0 which gives terms of the same order in the upper bound on the error of approximation. For a given initial condition X 0 =(W 0 , λ 0 , Ξ 0 )∈D 1 and a time of escape T 1 > 0, we assume that the solution X . (t) = X . i (t) i∈{1,...,8} = (W . (t), λ (t), Ξ . (t)) governed by Ξ + H , does not escape of D 1 for all |t| ≤ T 1 .