This is a preprint of: “On co-orbital quasi-periodic motion in the three-body problem”, Josep Maria Cors, Jesús F. Palacián, Patricia Yanguas, SIAM J. Appl. Dyn. Syst., vol. 18(1), 334–353, 2019. DOI: [10.1137/18M1190859] On co-orbital quasi-periodic motion in the three-body problem Josep M. Cors∗ Jesús F. Palacián† Patricia Yanguas‡ Abstract Within the framework of the planar three-body problem we establish the existence of quasi- periodic motions and KAM 4-tori related to the co-orbital motion of two small moons about a large planet where the moons move in nearly circular orbits with almost equal radii. The approach is based on a combination of normal form and symplectic reduction theories and the application of a KAM theorem for high-order degenerate systems. To accomplish our results we need to expand the Hamiltonian of the three-body problem as a perturbation of two uncoupled Kepler problems. This approximation is valid in the region of phase space where co-orbital solutions occur. MSC Classifications: 70F10, 70K65, 37J40, 37N05 Key Words: Three-body problem; symplectic scaling; co-orbital regime; 1:1 mean-motion resonance; normalization and reduction; KAM theory for multiscale systems; quasi-periodic motion and invariant 4-tori. 1 Introduction The purpose of the paper is to prove the existence of quasi-periodic motions related to the co-orbital motion of two small moons orbiting a big planet by means of KAM theory. We refer to co-orbital motion when two bodies are in a 1:1 mean-motion resonance with respect to each other while orbiting a central body. A horseshoe orbit is a type of co-orbital motion where two bodies in a 1:1 mean-motion resonance maintain the same average distance from the central body in libration around points 180◦ apart. The horseshoe motion and the co-orbital dynamics have been studied both analytically and numerically. One approach deals with the restricted three-body problem, assuming that one of the small bodies has negligible mass, see for instance [9, 2]. On the other hand, co-orbital motion can be approximated by two independent solutions of a two-body ∗ Departament de Matemàtiques, Universitat Politècnica de Catalunya, cors@epsem.upc.edu † Departamento de Estadı́stica, Informática y Matemáticas and Institute for Advanced Materials (INAMAT), Universidad Pública de Navarra, 31006 Pamplona, Spain, palacian@unavarra.es ‡ Departamento de Estadı́stica, Informática y Matemáticas and Institute for Advanced Materials (INAMAT), Universidad Pública de Navarra, 31006 Pamplona, Spain, yanguas@unavarra.es 1 problem, one when the small bodies do not interact, and another when the interaction between the small bodies dominates the motion. Finally, a matching of both solutions completes the process, as it is done in [20]. In the framework of the three-body problem, usually, the interaction between the two small bodies is treated as a perturbation of decoupled Kepler problems, see for example, [11, 12]. Even though, when the co-orbital motion is considered one has to keep in mind that the distance between the orbits of the small bodies is also small and their interaction will become large when the two small bodies are close. Due that, there are three small quantities to consider. Following the ideas in [5], we assume that the difference between the semi-major axes of the orbits when the small bodies are far apart is of order ε compared with the average radius of their orbits. Second, we assume that the masses of the small bodies are of order εa compared to the big mass or central body. Finally, since we are looking for orbits that has no syzygy when the small bodies are in the same side of the central body, we asume that the minimum angular separation, if we look at the particles from the big mass, is of order εb . In order to prepare the three-body problem so that we can apply an adequate KAM theorem we need to perform a sequence of symplectic transformations that lead us to the co-orbital regime. The transformations involving the small parameter are in the setting of the symplectic scaling techniques of Meyer, that he has applied successfully to the N -body problem, with the aim of accounting for different regimes of the problem, see [12]. These techniques were already applied in [5], and the authors obtained some relationships between the parameters involved in the co-orbital motion. The constraints concerning the parameters are crucial in our approach as they allow us to arrange the Hamiltonian function properly, placing its terms at different orders of the perturbation in an accurate way. Besides, we have applied many of the transformations of [5], some of them slightly generalized, to perform our procedure. A normal form transformation of the system is carried out assuming that the zeroth order of the Hamiltonian is in 1:1 resonance. After truncating higher order terms we apply symplectic reduction theory and analyze the resulting reduced system. The theory of Hamiltonian reduction goes back to Laplace and Lagrange although in its modern version it is usually attributed to the works of Reeb [18], Moser [17], Meyer [11] and Marsden and Weinstein [10] and is known nowadays as Meyer-Marsden-Weinstein reduction or regular reduction. However, often the reduced space is not a manifold but presents singularities, thus the theory lies in the context of the so called singular reduction [6]. This is indeed the case of our research. The zone in phase space where the co-orbital motion occurs is discussed conveniently and some coordinates are introduced so that the Hamiltonian is organized adequately and expressed in action- angle coordinates. As the action variables we have built arise in the Hamiltonian function at three different scales, showing the highly degenerate character of the Hamiltonian, usual KAM theorems, even those dealing with situations where the degeneration is removed, cannot be applied. Thus we apply a result due to Han, Li and Yi [8] that is suitable for our degenerate system. As far as we know there are not similar studies dealing with the persistence of invariant tori for the co-orbital motion of the three-body problem. In [19] the authors apply averaging theory associated to the 1:1 resonance in the context of the planar three-body problem, also focusing on the co-orbital regime. All their results are only related to the system obtained after truncating higher order terms. The paper is organized as follows. In Section 2 we provide the equations of motion of the three- body problem as well as the scalings and successive symplectic changes that are needed to focus our attention on the co-orbital motion. The purpose of Section 3 is to normalize the system with 2 respect to an angle that accounts for the 1:1 resonance. This allows us to simplify drastically the Hamiltonian function and prepare it, after truncation of the remainder, to apply reduction theory. In Section 4 the two degrees of freedom reduced system is analyzed, identifying essentially three types of motion, namely the ones related to the Lagrangian and Eulerian solutions of the three- body problem plus the motions of co-orbital nature that are the main concern of our approach. KAM theory is applied in Section 5 using the theorem of Han, Li, Yi [8] and applying some more transformations. It leads to our main result, that is Theorem 5.1. Finally some conclusions and prospective work are commented in Section 6. 2 Equations and scalings In this section, we derive the equations of motion for two small co-orbital particles in near coplanar circular orbits around a big mass, as a perturbation of two uncoupled Kepler problems. For bet- ter understanding and self-contained reading we detail the most relevant symplectic co-ordinates transformations given in [5]. We start by considering the Hamiltonian of the three-body problem p20 p2 p2 m1 m2 m1 m2 H(q0 , q1 , q2 , p0 , p1 , p2 ) = + 1 + 2 − − − (1) 2 2m1 2m2 |q1 − q0 | |q2 − q0 | |q2 − q1 | where q0 , p0 are the position and momentum in the plane of a big mass M = 1, qi , pi , i = 1, 2 are the positions and momenta of the small particles with masses mi , i = 1, 2, respectively, where the distance units have been chosen such that the gravitational constant G = 1. First we perform a symplectic change of variables (q, p) → (u, v), like Jacobi variables [13], where the two small particles are treated equivalently and u3 , v3 represent the center of mass and total momentum, respectively. Clearly, the change of variables isolates the total momentum and the center of mass, and so, it suffices to consider only variables u1 , u2 , v1 , v2 . Let us introduce a small parameter related to the masses, through mi = εa µi , i = 1, 2, where µ1 , µ2 are of the same order of the big mass M = 1. Simultaneously, we rescale the momentum variables (v1 , v2 ) with the same factor vi = εa wi , i = 1, 2. Notice that the change of variables from (u, v) → (u, w) is symplectic with multiplier ε−a . Next we perform the usual change to rotating polar coordinates (r, θ, R, Θ). The new Hamiltonian is     1 2 Θ21 1 2 Θ22 µ1 µ2 H(r, θ, R, Θ) = R1 + 2 + R2 + 2 − Θ1 − Θ2 − − 2µ1 r1 2µ2 r2 r1 r2 2 2 2 −µ1 µ2 R1 R2 Θ1 Θ22 +εa p 2 + + + + (2) r − 2r 1 r2 cos(θ2 − θ1 ) + r 2 2 2 2r12 2r22  1  2   Θ1 Θ2 R2 Θ1 R1 Θ2 + + R1 R2 cos(θ2 − θ1 ) + − sin(θ2 − θ1 ) . r1 r2 r1 r2 At this point we introduce the second small parameter, the difference between the semimajor axes. We assume that this difference is of order ε, and perform a symplectic change of variables with multiplier ε−1 given by ερ1 = r1 − 1, εΦ1 = Θ1 − µ1 , (3) ερ2 = r2 − 1, εΦ2 = Θ2 − µ2 . 3 We stress that in order to avoid the singularity that appears at ε = 0 we rescale the time appropriately. We expand and collect terms with respect to powers of ε, with the goal of getting a simple system at order zero that reflects the co-orbital motion. Due that, on the one hand, we have to keep terms of order ε2 and lower to represent accurately the motion of the small bodies when they are far apart. Moreover, we have to select upper bounds for a and conditions on b, where θ, 2π − θ ≥ O(εb ), such that finitely many terms, resulting from the interaction between the small bodies when they are close to one another, are of order less than or equal to 2 in ε, as well. This is precisely the main result stated in [5] where the relative size of the small quantities such that the planar three-body problem admits co-orbital motion is given. Theorem 2.1 (Cors-Hall). Let ε be the ratio of the difference between the radii of the orbits of the small moons with the average radius, and suppose the ratio of the moons’ masses to the mass of the planet is of order εa and the minimum angular separation of the moons is of order εb . If 2 < a < 5/2, then the system admits co-orbital motion, and in this case b = a − 2. Following the previous theorem, when a ∈ (2, 52 ), expanding and collecting similar power terms of ε, we obtain the following Hamiltonian expressed in variables (ρ, θ, R, Φ): H = H0 + ε2 H2 + εa Ha + ε3 H3 + εa+1 Ha+1 + ε4 H4 + O(εa+2 ) (4) where R12 R2 H0 = + 2, 2µ1 2µ2   1 Φ21 Φ22 2 2 H2 = + + µ1 ρ1 + µ2 ρ2 − 4(Φ1 ρ1 + Φ2 ρ2 ) , 2 µ1 µ2 R12 R22 Ha = + + (µ1 R2 − µ2 R1 ) sin(θ2 − θ1 ) + (µ1 µ2 + R1 R2 ) cos(θ2 − θ1 ) 2 2   µ1 µ2 θ2 − θ1 + csc , 2 2 Φ21 ρ1 Φ22 ρ2 H3 = − − − µ1 ρ31 − µ2 ρ32 + 3(Φ2 ρ22 + Φ1 ρ21 ), µ1 µ2   2 2 1 θ2 − θ1 Ha+1 = − µ1 ρ1 − µ2 ρ2 + µ1 Φ1 + µ2 Φ2 − µ1 µ2 (ρ1 + ρ2 ) csc 4 2 + (µ2 R1 ρ2 − µ1 R2 ρ1 − Φ2 R1 + Φ1 R2 ) sin(θ2 − θ1 ) + (µ2 Φ1 + µ1 Φ2 − µ1 µ2 (ρ1 + ρ2 )) cos(θ2 − θ1 ),   3 Φ21 ρ21 Φ22 ρ22 H4 = + + µ1 ρ1 + µ2 ρ2 − 4(Φ2 ρ32 + Φ1 ρ31 ). 4 4 2 µ1 µ2 We remark that the order of the power terms in ε obtained does not depend on the admissible values of a. The most important term resulting from the interaction of the two small bodies is given by csc((θ2 − θ1 )/2) in Ha . Clearly, at the moment of minimum angular separation, that is, when θ or 2π − θ are equal to εb , it will be of order ε2 . Finally, taking into account the forthcoming changes of variables we have to keep here at least terms to order ε4 . 4 Since the dependence of the Hamiltonian with respect to the two angles θ1 and θ2 is only given by their difference θ2 − θ1 we perform a symplectic transformation with multiplier µ11 + µ12 that reduces the four degrees of freedom system onto a three degrees of freedom system, as follows q ρ = ρ2 − ρ1 , ρ̄ = µµ12 (ρ1 + µµ12 ρ2 ), q θ = θ2 − θ1 , θ̄ = µµ12 (θ1 + µµ21 θ2 ), (5) R = Rµ22 − Rµ11 , R̄ = R√1µ+R 2 1 µ2 , Φ2 Φ1 Φ √1 +Φ2 . Φ= µ2 − µ1 , Φ̄ = µ1 µ2 Finally, we separate the fast and slow variables, in order to emphasize the 1:1 resonance, by setting 1 1 ε− 2 γ = ρ, ε− 2 γ̄ = ρ̄, 1 1 (6) ε 2 P = R, ε 2 P̄ = R̄, and rescale time to obtain the Hamiltonian H(γ, γ̄, θ, θ̄, P, P̄ , Φ, Φ̄) to fifth order in  = ε1/4 and a = 49 : H = K0 + 2 K2 + 4 K4 + 5 K5 + O(6 ), (7) where µ21 µ22 2 2 2 2  K0 = P + P̄ + γ + γ̄ , 2(µ1 + µ2 )2 µ21 µ22 √ 2 2  3  K2 = − µ 1 µ 2 γ̄ 3γ + γ̄ + 2(µ 1 + µ 2 )(γΦ + γ̄ Φ̄) + (µ 1 − µ 2 )γ , (µ1 + µ2 )3 µ21 µ22   K4 = − 4 µ21 Φ2 + 6Φγ 2 + Φ̄2 + 3γ 4 + µ22 Φ2 − 6Φγ 2 + Φ̄2 + 3γ 4 2(µ1 + µ2 ) 3/2 √   √ 3/2   +6µ1 µ2 2γγ̄ Φ + γ 2 + Φ̄ γ 2 + γ̄ 2 + 6 µ1 µ2 2γγ̄ Φ − γ 2 + Φ̄ γ 2 + γ̄ 2   +µ1 µ2 2Φ2 + 2Φ̄2 − 3γ 4 + 18γ 2 γ̄ 2 + 3γ̄ 4 − (µ1 + µ2 )3 csc(θ/2) , µ21 µ22 K5 = cos θ. µ1 + µ2 Although the factors that appear in the different orders in  depend on the admissible chosen value of a, this will not affect future conclusions. As a compromise we have taken a = 9/4, the midpoint of the admissible interval. On the other hand, K2 contains terms coming from H2 and H3 . The terms of K4 proceed from H2 , Ha , H3 and H4 . Finally, K5 only collects terms from Ha . 3 Normalization To achieve the normalization process it is convenient to introduce standard action-angle variables, namely J, ψ, p p γ = 2(J1 − J2 ) sin ψ1 , γ̄ = 2(J1 − J2 ) cos ψ1 , √ √ P = 2J2 sin(ψ1 + ψ2 ), P̄ = 2J2 cos(ψ1 + ψ2 ), where J1 ≥ J2 ≥ 0. 5 Expressing (7) in these variables, the zeroth-order term K0 after a suitable rescaling of time is simply J1 , while the rest of terms are transformed and scaled accordingly. Indeed the perturbative terms become a finite Fourier series in ψj whose coefficients are functions of Jk , that also depend on Φ and Φ̄. In addition to that one has to consider the terms related to θ in K4 and K5 that are not altered by this last transformation. Then, we average the resulting Hamiltonian with respect to the angle ψ1 in order to obtain a two degrees of freedom Hamiltonian, after truncating higher order terms. The averaging process is customary and is based on the Lie transformations approach for Hamiltonian systems, see [7]. We have used the software Mathematica to carry out the computations since the process involves the calculation of many terms. Alternatively one could have performed the normalization process in complex-rectangular coor- dinates derived from γ, γ̄, P , P̄ , taking into account that K0 is in 1:1 resonance. Both procedures yield the same result. After some simplifications we end up with the normalized Hamilton function given by 34  H = J1 − (µ21 + µ22 )(J1 − J2 )2 + (µ1 + µ2 )2 (Φ2 + Φ̄2 ) 2(µ1 + µ2 )2 √ +2(µ21 − µ22 )(J1 − J2 )Φ − µ1 µ2 (J12 − 6J1 J2 + 4J22 ) + 2 µ1 µ2 (µ1 + µ2 )J1 Φ̄ √ p  (8) +4 µ1 µ2 (J1 − J2 )J2 (µ1 − µ2 )(J1 − J2 ) + (µ1 + µ2 )Φ cos ψ2  −6µ1 µ2 J2 (J1 − J2 ) cos(2ψ2 ) − (µ1 + µ2 )3 csc(θ/2) +5 (µ1 + µ2 ) cos θ + O(8 ). We stress that the terms of orders 6 and 7 are identically zero and the next significant term occurs at order 8 in , that is, at order 2 in ε. The explicit expression of O(8 ) has also been computed as it is needed in the determination of the stability character of the triangular and collinear solutions studied in Section 4. An important remark is that the normalization we have carried out induces the introduction of a formal first integral, namely J1 , because we are averaging with respect to a single angle, thus it does not involve the management of small divisors. In turn it implies that the terms of the remainder can be considered exponentially small with respect to the small parameter , revealing a good symptom of the quality of the normal form approach. For more details, see for instance Chapter 6 of [1] and the recent paper [3]. We also notice that if the remainder is dropped, the equations of motion associated to H in (8) define a Hamiltonian system of two degrees of freedom, given by the action-angle pairs ψ2 /J2 and θ/Φ. 4 Reduction and analysis of the reduced system 4.1 Reduced space After truncating the terms of order eight and higher in , the normal form Hamiltonian H can be reduced with respect to the symmetry introduced through the normalisation process. Besides, it is also reduced by the symmetry induced by the integral Φ̄. 6 To achieve this step we define suitable global variables for the reduced two degrees of freedom system. These coordinates are denoted by σi and ai , i = 1, . . . , 4 and they are generically called invariants. Regarding the approximate integral of motion introduced through the normalisation process, that is J1 , it is clear that since it is related to the interaction occurring between the pairs J1 /ψ1 and J2 /ψ2 , the σi are the ones corresponding to the 1:1 resonance. Thence they are introduced in terms of γ, γ̄, P and P̄ through σ1 = γ 2 + P 2 , σ2 = γ̄ 2 + P̄ 2 , σ3 = γγ̄ + P P̄ , σ4 = γ P̄ − γ̄P, (9) see for instance [6] or [16]. The invariants associated to the reduction due to the integral Φ̄ are obtained following the procedure given in [14, 15]. Since this integral is essentially the one dealing with the conservation of the angular momentum vector there are rectangular symplectic coordinates, say x1 , x2 , y1 , y2 , such that Φ̄ is identified with the Coriolis term x2 y1 − x1 y2 and the corresponding invariants related to it become a1 = x21 + x22 , a2 = y12 + y22 , a3 = x1 y1 + x2 y2 , a4 = x2 y1 − x1 y2 . Then, the transformation x1 = √1 (u1 − v2 ), x2 = √1 (u2 − v1 ), y1 = √1 (u2 + v1 ), y2 = √1 (u1 + v2 ), 2 2 2 2 converts x2 y1 − x1 y2 into 1 2 (u22 + v22 ) − 21 (u21 + v12 ), putting in emphasis the fact that the Coriolis term can be interpreted as a harmonic oscillator in 1:−1 resonance. Next we introduce standard Poincaré action-angle coordinates Υi /ξi , i = 1, 2: p p p p u1 = 2Υ1 sin ξ1 , u2 = 2Υ2 sin ξ2 , v1 = 2Υ1 cos ξ1 , v2 = 2Υ2 cos ξ2 , and perform the linear transformation Υ1 = Φ, Υ2 = Φ + Φ̄, ξ1 = θ − θ̄, ξ2 = θ̄, ending up with p a1 = 2Φ + Φ̄ − 2 Φ(Φ + Φ̄) sin θ, p a2 = 2Φ + Φ̄ + 2 Φ(Φ + Φ̄) sin θ, p (10) a3 = −2 Φ(Φ + Φ̄) cos θ, a4 = Φ̄, for Φ ≥ 0 and Φ̄ ≥ 0 or for Φ ≥ −Φ̄ ≥ 0. Similar linear transformations can be defined for other values of Φ and Φ̄. Concretely, for Φ ≥ 0 and Φ̄ ≤ 0 or for Φ ≥ Φ̄ ≥ 0 we get p p a1 = 2Φ − Φ̄ − 2 Φ(Φ − Φ̄) sin θ, a2 = 2Φ − Φ̄ + 2 Φ(Φ − Φ̄) sin θ, p (11) a3 = −2 Φ(Φ − Φ̄) cos θ, a4 = −Φ̄, while when Φ ≤ 0 and Φ̄ ≥ 0 or when Φ ≤ Φ̄ ≤ 0 we define p p a1 = −2Φ + Φ̄ + 2 Φ(Φ − Φ̄) sin θ, a2 = −2Φ + Φ̄ − 2 Φ(Φ − Φ̄) sin θ, p (12) a3 = −2 Φ(Φ − Φ̄) cos θ, a4 = Φ̄, 7 and for Φ ≤ 0 and Φ̄ ≤ 0 or for Φ ≤ Φ̄ ≤ 0 we apply p p a1 = −2Φ − Φ̄ + 2 Φ(Φ + Φ̄) sin θ, a2 = −2Φ − Φ̄ − 2 Φ(Φ + Φ̄) sin θ, p (13) a3 = −2 Φ(Φ + Φ̄) cos θ, a4 = −Φ̄. Fixing the integrals J1 and Φ̄ by means of J1 = j1 > 0 and Φ̄ = C ∈ R we arrive at the following constraints, two of them involving the σi and the other two for the aj . Specifically, we get σ1 + σ2 = 2j1 , σ1 (2j1 − σ1 ) = σ32 + σ42 , a23 + C 2 = a1 a2 , a4 = ±C. (14) The equations (14) define the reduced space Bj1 ,C . It is readily deduced that Bj1 ,C is the Cartesian product of a two-dimensional sphere and one sheet of a two-sheet hyperboloid of revolution. A picture is given in Fig. 1. × Figure 1: Reduced space Bj1 ,C . The reduced Hamiltonian is related to the truncation at order 5 of H defined in (8). It is given in terms of the invariants and we call it H̄(σi , ai ; j1 , C). The phase space associated to H̄ is Bj1 ,C . The two-sheet hyperboloid is a regular surface provided C 6= 0. More specifically, the closest point of the hyperboloid to the origin of the a-space is reached for Φ = 0 and has coordinates a1 = a2 = a4 = C, a3 = 0. This point becomes singular when C = 0, thus in this case the reduction process lies in the setting of singular reduction, see [6] and also [14, 16]. Nevertheless in our particular setting when C = 0 the relevant points are the regular ones in Bj1 ,0 . Therefore the results (mainly of local nature) of regular reduction apply, see for instance [21]. For convenience we will work with the coordinates Φ, θ instead of the invariants ai . The reason is that the forthcoming expressions become shorter using the action-angle pair and provide a better geometrical interpretation of the motions we are studying. Thus, we notice that θ is well defined as soon as Φ 6= 0. Moreover we need to exclude the values θ = 0, 2π due to the cosecant terms in the perturbation, but this is compatible with the fact that θ ∈ [εb , 2π − εb ] where b = a − 2 = 1/4 for the choice of a made above. 8 4.2 Relative equilibria The equations of motion are built as σ̇i = {σi , H̄}, ȧi = {ai , H̄}. (15) To compute these equations in an explicit way one needs to use the Poisson parentheses {σ1 , σ2 } = 0, {σ1 , σ3 } = 2σ4 , {σ1 , σ4 } = −2σ3 , {σ2 , σ3 } = −2σ4 , {σ2 , σ4 } = 2σ3 , {σ3 , σ4 } = σ1 − σ2 . Analogously, the parentheses involving the ai are {a1 , a2 } = 4a3 , {a1 , a3 } = 2a1 , {a2 , a3 } = −2a2 and the rest are zero; see also the details in [14]. Besides, {σi , aj } = 0 also holds. We have performed the complete analysis of the existing relative equilibria corresponding to (15), discarding the ones related to the singular point (ai = 0, C = 0) in the hyperboloid as this point has nothing to do with the co-orbital regime. In our analysis we have taken into account the constraints (14), ending up with the existence of six different critical points. However, instead of presenting explicitly the coordinates of these points in the invariants, we have preferred to pass to the action-angle coordinates ψ2 /J2 , θ/Φ. For this we use the fact that J1 = 21 (σ1 + σ2 ), J2 = 21 σ2 , together with σ4 σ3 sin ψ2 = − , cos ψ2 = , σ32 + σ42 σ32 + σ42 while the expressions of θ, Φ are readily deduced from (10), leading to   −1 a1 − a2 θ = ± tan , 2a3 1 1 Φ = (a1 + a2 ∓ C) when Φ ≥ 0 and Φ = − (a1 + a2 ∓ C) when Φ ≤ 0. 4 4 The equilibria (ψ2i , θi , J2i , Φi ), i = 1, . . . , 6 are given by       π µ2 j1 µ1 j1 µ2 j1 µ1 j1 5π µ2 j1 µ1 j1 0, , ,− , 0, π, ,− , 0, , ,− , 3 µ1 + µ2 µ1 + µ2 µ1 + µ2 µ1 + µ2 3 µ1 + µ2 µ1 + µ2       π µ1 j1 µ2 j1 µ1 j1 µ2 j1 5π µ1 j1 µ2 j1 π, , , , π, π, , , π, , , , 3 µ1 + µ2 µ1 + µ2 µ1 + µ2 µ1 + µ2 3 µ1 + µ2 µ1 + µ2 where in the first row we have written sequentially (ψ2i , θi , J2i , Φi ), i = 1, 2, 3 and in the second row, (ψ2i , θi , J2i , Φi ), i = 4, 5, 6. These equilibria are related to the Eulerian (collinear) and Lagrangian (triangular) solutions of the planar three-body problem. In particular, the two points with θi = π correspond to the collinear motions whereas the rest of the points are the ones related to the Lagrange triangular solutions. Indeed, for the collinear solution (ψ22 , θ2 , J22 , Φ2 ) the configuration is m1 -M -m2 while the point (ψ25 , θ5 , J25 , Φ5 ) has configuration m2 -M -m1 . We remark that the big body has to be always in the middle of the small bodies (equivalently θ cannot be zero or 2π) as otherwise it would lead to an inadmissible configuration. For the Lagrangian points, the bodies with masses M , m1 , m2 form an equilateral triangle. More precisely the point (ψ21 , θ1 , J21 , Φ1 ) represents the configuration with the body of mass M at 9 the bottom of the triangle, whereas the body of mass m1 is placed on the left vertex and the body of mass m2 on the right vertex. The point (ψ24 , θ4 , J24 , Φ4 ) has the M mass at the bottom of the triangle while the body of mass m2 is on the left vertex and the one whose mass is m1 on the right one. Besides, the point (ψ23 , θ3 , J23 , Φ3 ) represents the configuration with the body of mass M is at the top of the triangle, the body whose mass is m1 appears on the left vertex and the one with mass m2 on the right one. Finally the point (ψ26 , θ6 , J26 , Φ6 ) has the body of mass M placed at the bottom, the body of mass m2 is on the left vertex and the one with mass m1 on the right vertex. Trying to obtain the stability character of the six critical points, we notice that all lead to degenerate situations. That is, in the six cases the 4 × 4-matrices associated to the linearizations have a pair of null eigenvalues, thus the Hessians with respect to the coordinates (ψ2 , θ, J2 , Φ) evaluated at the equilibria are zero. To resolve this degeneracy one has to push the calculations to higher order, in particular normalizing explicitly the terms of order 8 . For this, having in mind the choice of a, one has to include terms of orders a + 2 and 5 in ε of the Hamiltonian (4). Thus, the successive changes of coordinates detailed in Sections 2 and 3, including the normalization, are carried out up to order 8 in . The normalized Hamiltonian is too big to be reproduced here and so, we do not write it down ex- plicitly. However, the coordinates for the equilibria do not vary with respect to the 5 -approximation. The eigenvalues of the linearization associated with (ψ2i , θi , J2i , Φi ), i = 1, 3, 4, 6 after expanding them in powers of  up to order 8 is √ 75ı8 µ2 j12 3 3ı 9/2 √ ± , ±  µ1 + µ2 , 4(µ1 + µ2 )2 2 concluding parametric stability for the Lagrangian points. These equilibria are of the type centre × centre. This result should be expected as the masses m1 , m2 are very small compared to the big body, and this situation corresponds to the case of linearly stable Lagrangian configurations in the three-body problem. Similarly, the eigenvalues of the linearization related to (ψ2i , θi , J2i , Φi ), i = 2, 5, after expansion to order 8 are √ 75ı8 µ2 j12 3 21 9/2 √ ± , ± √  µ1 + µ2 , 4(µ1 + µ2 )2 2 2 thus leading to instability for the collinear ones, as they are points of the form centre × saddle. Taking into consideration the 8 -approximation the dynamics corresponding to the six relative equilibria can be reconstructed to the unreduced system, that is, to the full Hamiltonian with four degrees of freedom. This is possible because the reduction process involves an approximate symmetry and an exact one, thus we can apply Reeb’s theory [18, 21] because we have made the normal form transformation averaging with respect to a single angle. Hence, we conclude that (ψ2i , θi , J2i , Φi ), i = 1, 3, 4, 6 correspond to invariant elliptic 2-tori of the full Hamiltonian system in four degrees of freedom. Besides, the equilibria (ψ2i , θi , J2i , Φi ), i = 2, 5 lead to invariant elliptic- hyperbolic 2-tori of the full Hamiltonian system in four degrees of freedom. 4.3 Dynamics corresponding to quasi-periodic co-orbital motions In order to identify in the reduced system the regime that is related to the quasi-periodic motions explaining the co-orbital motion between the small moons, we have to restrict ourselves to a specific 10 part of the reduced space corresponding to the pair θ/Φ, that is, to a piece of the two-sheet hyperboloid of revolution where co-orbital motion can be detected, “freezing” the dynamics on the σ-sphere. As the time variation of Φ is of order O(4 ), this action evolves slowly compared to the other actions, thus it can be considered as a slow variable. In a first approximation we can think of it as almost a constant, i.e. an adiabatic invariant and concentrate on the dynamics on the sphere, that is looking for the “equilibria” on the σ-sphere. Specifically, considering the subsystem σ̇i = {σi , H} one gets two equilibria, namely  √  + + + + 2µ1 j1 2µ2 j1 2 µ1 µ2 j1 e+ = (σ1 , σ2 , σ3 , σ4 ) = , , ,0 , µ1 + µ2 µ1 + µ2 µ1 + µ2  √  − − − − 2µ2 j1 2µ1 j1 2 µ1 µ2 j1 e− = (σ1 , σ2 , σ3 , σ4 ) = , ,− ,0 . µ1 + µ2 µ1 + µ2 µ1 + µ2 Indeed, there are more equilibria but they are not isolated points. The coordinates J2 , ψ2 corre- sponding to e+ , e− are ψ2+ = 0, J2+ = µ2 j1 /(µ1 + µ2 ) and ψ2− = π, J2− = µ1 j1 /(µ1 + µ2 ). If Φ ≈ Φ0 (a constant), the eigenvalues associated to e+ and e− are, respectively, 3ı (µ1 j1 + (µ1 + µ2 )Φ0 ) 3ı (−µ2 j1 + (µ1 + µ2 )Φ0 ) ± , ± . µ1 + µ2 µ1 + µ2 Thus, provided µ1 j1 + (µ1 + µ2 )Φ0 does not vanish, e+ corresponds to a centre and when −µ2 j1 + (µ1 + µ2 )Φ0 is not null, e− is a centre as well. Now we replace in the normal form Hamiltonian (8) the values of ψ2 , J2 corresponding to e+ and e− and focus on the dynamics related to θ/Φ. After some simplifications the relevant terms of the Hamiltonians become 3j1 √ 3  he+ = − (µ1 j1 + 2 µ1 µ2 C + 2µ1 Φ) − (C 2 + Φ2 ) + (µ1 + µ2 ) (2 cos θ − csc(θ/2)) , 2(µ1 + µ2 ) 2 2 3j1 √ 3  he− = − (µ2 j1 + 2 µ1 µ2 C − 2µ2 Φ) − (C 2 + Φ2 ) + (µ1 + µ2 ) (2 cos θ − csc(θ/2)) . 2(µ1 + µ2 ) 2 2 The analysis of the one degree of freedom Hamiltonians he+ , he− in terms of θ, Φ, for all values of the constants, yields a saddle point to which a separatrix is attached that surrounds two centres. Being more precise, the coordinates of the centres for e+ are θ = π/3, 5π/3 while Φ = −µ1 j1 /(µ1 + µ2 ). For the saddle one has θ = π with the same Φ as before. In case of e− one arrives at the same values of θ as in e+ but Φ = µ2 j1 /(µ1 + µ2 ). Then e+ refers to motions with negative Φ whereas e− is related to motions with positive Φ. A typical picture corresponding to Hamiltonian he− is depicted in Fig. 2. From this study it becomes clear that the saddle points are related to the collinear solutions whereas the centres correspond to the Lagrangian motions. In both cases they have to be understood as the projections of this type of motions on the plane θ-Φ. From the previous paragraph and also from Fig. 2, it is readily deduced that the quasi-periodic motions associated to the co-orbital regime can be identified by the libration regime in the outer side of the separatrix. In order to have a notion on the values that Φ takes for co-orbital motions, we compute the maximum and minimum values of Φ on the separatrix. Notice that they are reached 11 0.7 0.6 Φ 0.5 0.4 0.3 0 1 2 3 4 5 6 θ Figure 2: Projection of the dynamics onto the plane θ-Φ related with e− . when θ = π/3, 5π/3. Also, the values of the energy on the separatrix are obtained by plugging in the expressions of he+ , he− the coordinates of θ, Φ corresponding to the saddle. Next, these values are equated with the energies he+ , he− where we have replaced θ by π/3, arriving at two equations of degree 2 in Φ. The roots are the desired values for Φ. If θ = π/3, 5π/3, we find that Φ has to be outside the interval " r r # µ1 j1 2 µ1 j1 2 − − (µ1 + µ2 ) , − + (µ1 + µ2 ) (16) µ1 + µ2 3 µ1 + µ2 3 for e+ , while for e− we find that Φ should not belong to the interval " r r # µ2 j1 2 µ2 j1 2 − (µ1 + µ2 ) , + (µ1 + µ2 ) . (17) µ1 + µ2 3 µ1 + µ2 3 Notice that in accordance with the intervals obtained above, Φ has to be relatively far from zero. These intervals give approximate bounds for Φ because they have been obtained from the Hamil- tonian normal form after plugging in it the coordinates of the “equilibria” e+ , e− . However, they are good approximations in the sense that in the co-orbital regime ψ2 , J2 are not fixed but they are near their values ψ2± , J2± . The bounds for Φ should be taken into account when dealing with the applications to specific examples. Furthermore they are relevant in order to estimate the size of this action when applying the KAM theorem of Section 5. To accomplish the analysis of the quasi-periodic motions we need coordinates associated to the pair ψ2 /J2 such that they account for the motion in the σ-sphere in a neighbourhood of the points e+ , e− . The pair θ/Φ is allowed to vary in the outer side of the separatrix, i.e. in the libration regime of the projection onto the plane θ-Φ. 12 5 KAM 4-tori In this section we prove the existence of invariant 4-tori associated to the co-orbital motion. So, we focus in the region of phase space that is outside the separatrix in Fig. 2. We start by preparing Hamiltonian (8) in order to apply a KAM-type theorem. In our case, as the system is highly degenerate, with the actions appearing for their first time at three different scales, we use a theorem by Han, Li and Yi [8]. In order to obtain the quasi-periodic motions we need four actions and four angles. We already have an action-angle pair in (8), that is Φ/θ. Taking into account that the last step in the reduction procedure has been fixing J1 = j1 we can undo the reduction and consider J1 as the second action in our tori, with its conjugate angle ψ1 . With the goal of getting the third action-angle pair, we construct local symplectic coordinates (x, y) around the points e+ and e− appearing in the σ-sphere of the previous section. These coordinates are related to ψ2 , J2 and are variables associated to the 1:1 resonance, see [16]. The expression of the invariants (σ1 , . . . , σ4 ) as functions of (x, y) is: σ1 = x2 + y 2 , σ2 = 2j1 − (x2 + y 2 ), p σ3 = −x 2j1 − (x2 + y 2 ), p σ4 = y 2j1 − (x2 + y 2 ). Then, the points e+ and e− are transformed into s ! s ! 0 0 2µ1 j1 2µ2 j1 (x1 , y1 ) = − ,0 and (x02 , y20 ) = ,0 . µ1 + µ2 µ1 + µ2 From here on most of the computations are presented for (x01 , y10 ) although the construction for (x02 , y20 ) would be similar, emphasizing that the choice (x01 , y10 ) implies that Φ < 0, whereas (x02 , y20 ) is related to positive Φ. We shift the origin to (x01 , y10 ) by setting x = x̄ + x01 , y = ȳ + y10 . Next, we restrict to a region in phase space where |Φ| is big enough to be able to uncouple the action Φ from the corresponding angle θ. For that, we introduce a new small parameter, η, by means of  = η 8 and replace Φ by Φ/η 3 . This assumption should be in agreement with the bounds provided by the interval given in (16) for θ = π/3, 5π/3. In other words, for e+ we know that Φ has to be taken in such a way that Φ∗ < Φ < 0 where Φ∗ is the right endpoint of the interval (16), but such that |Φ| is big enough so that the scaling made above makes sense. Similarly, for e− we have to consider 0 < Φ < Φ# , where Φ# is the left endpoint of the interval (17) and Φ should be big enough. Applying all these changes to Hamiltonian (8) and expanding the resulting function around x̄ = ȳ = 0 we get   3 26 2 29 3µ1 H = J1 − η Φ + η − J1 Φ + f1 (J1 , Φ, x̄, ȳ) + O(η 32 ) (18) 2 µ1 + µ2 and the first appearance of θ is at order 32 in η through the cosecant term. We do not give the explicit expression of f1 because it is lengthy and does not provide significant information. Since we 13 are considering motions in the σ-sphere close to e+ , the coordinates x̄, ȳ are small and we expand f1 with respect to them up to degree twelve. At this point we define an action K and an angle κ associated to (x̄, ȳ) by means of the canonical change of variables s r 2µ2 2(µ1 + µ2 ) x̄ = K sin κ, ȳ = K cos κ. µ1 + µ2 µ2 Applying this transformation to Hamiltonian (18) we get 3  µ1 H = J1 − η 26 Φ2 + 3η 29 − J1 Φ + ΦK 2 µ1 + µ2 µ1/2 ΦK 3/2  − µ1 + 2µ2 + µ1 cos(2κ) sin κ (19) 3/2 1/2 2µ2 J1  +f2 (J1 , Φ, K, κ) + O(η 32 ). Although it is calculated, we do not write down f2 explicitly, but it is a Fourier series in κ with coefficients depending on powers of K 1/2 starting at 2. With the goal of simplifying the Hamiltonian function (19) we normalize the terms of order O(η 32 ) with respect to the angle κ, arriving at   3 26 2 29 µ1 1 H = J1 − η Φ + 3η J1 (K − Φ) + (2Φ − K)K + O(η 32 ). (20) 2 µ1 + µ2 2 To achieve this normalization we have proceeded applying several steps of the Lie transformation and the dependence of the terms of order 32 on the angle κ is completely removed. Indeed because of the particular form of f2 it is enough to perform two steps of the normalizing transformation. Besides, proceeding in a similar fashion, for (x02 , y20 ) we get   3 26 2 29 µ2 1 H = J1 − η Φ + 3η J1 (Φ + K) − (2Φ + K)K + O(η 32 ). (21) 2 µ1 + µ2 2 At this point we have the Hamiltonian ready to apply Han-Li-Yi’s Theorem [8]. Theorem 5.1. For a fixed ε > 0, there are KAM 4-tori related with the co-orbital motions of the two small masses in the planar three-body problem (1) and the excluding measure for the existence of these invariant tori is of order O(εδ/32 ) with 0 < δ < 1/5. Proof. First we apply Han-Li-Yi’s Theorem [8] to Hamiltonian (20) to obtain the 3D tori. The application to (21) is analogous and leads to similar results valid for Φ > 0. The three actions we need to achieve the persistence of tori are J1 , Φ and K and the first terms of the Hamiltonian H are h0 (J1 ) = J1 , 3 h1 (J1 , Φ) = − Φ2 , (22) 2   µ1 1 h2 (J1 , Φ, K) = 3 J1 (K − Φ) + (2Φ − K)K . µ1 + µ2 2 14 We consider I = (J1 , Φ, K) and the frequency vector   ∂h0 ∂h1 ∂h2 Ω ≡ (Ω1 , Ω2 , Ω3 ) = , , , ∂J1 ∂Φ ∂K from where we compute the matrix  ∂Ω1 ∂Ω1 ∂Ω    1 Ω1 1 0 0 0  ∂J1 ∂Φ ∂K        ∂Ω2 ∂Ω2 ∂Ω2   −3Φ 0 −3 0 ,  Ω2 =     ∂J1 ∂Φ ∂K      J1 µ1 3µ1 ∂Ω3 ∂Ω3 ∂Ω3 3 Φ−K + 3 −3 Ω3 µ1 + µ2 µ1 + µ2 ∂J1 ∂Φ ∂K and obtain n o α Rank ∂I Ω(I) : 0 ≤ |α| ≤ N = 3, with N = 1. In the notation of Han-Li-Yi’s theorem, the orders where the actions appear are m1 = 26, m2 = 29; the number of actions hi in (22) depend on are n0 = 1, n1 = 2 and n2 = 3. Then, b = m1 (n1 − n0 ) + m2 (n2 − n1 ) = 55. Therefore, η b = ε55/32 , while the angles appear at order η 32 = ε (thus, at order 1 in ε) and as b = 55/32 > 1, the excluding measure is of the order εδ/32 with 0 < δ < 1/5. The same conclusion is achieved for Hamiltonian (21). The last step is undoing the reduction by the rotational symmetry by taking into account that C = Φ̄ is an action and has conjugate angle θ̄. Since the reduction related with Φ̄ is exact, that is, it does not require any truncation step, its reconstruction is immediate. In this way we obtain the KAM 4-tori. The corresponding quasi-periodic motions are of prograde nature since in (3) we have restricted ourselves to positive angular momentum Θ1 + Θ2 . However the construction leading to the persis- tence of invariant tori of retrograde type is essentially the same as the one exposed here. Remark 5.2. The estimate O(εδ/32 ) can be refined in order to get a much better measure of the tori that remain after adding the perturbation. This has to be achieved by pushing the occurrence of the angles θ and κ up to order 2 in ε, i.e. up to order 64 in η. For that, we first need to average Hamiltonian (18) with respect to θ. This process is delicate as the cosecant term in θ leads to logarithmic terms in the corresponding generating function. However, both the averaged Hamiltonian and the normalizing transformation are analytic functions in the allowed domain for θ. The second step is the averaging of the resulting Hamiltonian with respect to κ which is rather customary. However, a detail regarding this second normalization when compared with the one performed above, that is, the one leading to Hamiltonian (20), is that we have scaled the coordinates x̄, ȳ (and consequently the whole Hamiltonian) getting a slight improvement of the estimates. Finally, after making the double averaging we arrive at a Hamiltonian such that the dependence of the angle terms starts at order ε2 . Specifically for e+ we get 3 µ1 H = J1 − η 26 Φ2 − 3η 29 J1 Φ + 3η 31 ΦK + O(η 32 ), (23) 2 µ1 + µ2 and an analogous expression for e− . 15 Although in Hamiltonian (23) there are terms between the orders 32 and 64 in η they depend only on the actions. One has to take into account the form of the normalized Hamiltonian (8), where the terms of orders 6 and 7 in  were empty. Also, after the normalization over θ and κ performed in (23), the angles κ and θ appear for the first time at O(η 64 ), while ψ2 appears at order 72 in η. In this way, since 57/32 < 2, applying Remark (2) in [8] the final estimate of the excluding measure of the tori improves substantially becoming O(ε57/32 ). 6 Conclusions and future work We have developed a theory with the aim of accomplishing the persistence of quasi-periodic motions and related invariant 4-tori associated to the co-orbital motion of two small moons around a planet, in the setting of the planar three-body problem. An important point in our achievement is the introduction of a small parameter, taken as the ratio of the difference between the radii of the orbits of the small moons with the average radius. The choice of the small parameter together with the execution of successive changes of coordinates allow us to put the Hamiltonian function in an adequate form in order to apply averaging theory and reduce after truncating the related higher-order terms. For the reduced system we easily identify in the corresponding phase space the regime where co-orbital motion occurs. Since the perturbation of the Hamiltonian function comes out at different orders, a crucial step is the application of an appropriate theorem to carry out the persistence of the KAM tori in multiscale Hamiltonians, with at least three different scales. Indeed, standard theorems dealing with degenerate situations cannot be applied in our context. Moreover, a refinement in the normalization process, removing the dependence of the normal form Hamiltonian with respect to θ up to order O(ε2 ) has allowed us to improve the estimate on the excluding measure of the tori up to an order O(ε1.781... ). The present approach can be generalized by dropping the restriction of considering moderate to big angular momentum Φ. More specifically, one can define an adequate pair of action-angle coor- dinates by means of some elliptic functions, performing a careful analysis of the related expressions. In this context one can apply our analysis to the system formed by Saturn and its moons Janus and Epimetheus that move on a co-orbital regime. This is an ongoing study and will appear elsewhere. Finally, according to the data on Saturn’s moons provided by the Cassini mission [4], the orbits of Janus and Epimetheus are three-dimensional, thence it would make sense to perform an analysis of their motion in the framework of the spatial three-body problem. This study would need more reductions, enlarging therefore the present approach. Acknowledgments: The authors want to thank the Centre de Recerca Matemàtica (CRM) in Bellaterra, Barcelona, where this study was initiated in the framework of the Research Program on Central Configurations, Periodic Orbits and Beyond in Celestial Mechanics hosted at CRM from January to July, 2014. J. M. Cors was partially supported by grants MTM2016-77278-P (FEDER) and AGAUR grant 2017 SGR 1617. J. F. Palacián and P. Yanguas have been partially supported by grants MTM 2014-59433-C2-1-P and MTM 2017-88137-C2-1-P. 16 References [1] Arnold, V. I., Kozlov, V. V. and Neishtadt, A. I., Dynamical Systems III, Third Edition, Encyclopaedia of Mathematical Sciences 3, Springer, Berlin, 2006. [2] Barrabés, E. and Ollé, M., Invariant manifolds of L3 and horseshoe motion in the restricted three-body problem, Nonlinearity, 19 (9), 2006, 2065–2089. [3] Chartier, P., Murua, A. and Sanz-Serna, J. M., Higher-order averaging, formal series and numerical integration III: error bounds, Found. Comput. Math., 15 (2), 2015, 591–612. [4] Cooper, N. J., Renner, S., Murray, C. D. and Evans, M. W., Saturn’s inner satellites: orbits, masses, and the chaotic motion of Atlas from new Cassini imaging observations, Astron. J., 149 (27), 2015, 18 pp. [5] Cors, J. M. and Hall, G. R., Coorbital periodic orbits in the three body problem, SIAM J. Appl. Dyn. Syst., 2 (2), 2003, 219–237. [6] Cusham, R. H. and Bates, L. M., Global Aspects of Classical Integrable Systems, Second Edition, Birkhäuser, Basel, 2015. [7] Deprit, A., Canonical transformations depending on a small parameter, Celestial Mech., 1 (1), 1969, 12–30. [8] Han, Y., Li, Y. and Yi, Y., Invariant tori in Hamiltonian systems with high order proper degeneracy, Ann. Henri Poincaré, 10 (8), 2010, 1419–1436. [9] Llibre, J. and Ollé, M., The motion of Saturn coorbital satellites in the restricted three-body problem, Astronom. Astrophys., 378 (3), 2001, 1087–1099. [10] Marsden, J. and Weinstein, A., Reduction of symplectic manifolds with symmetry, Rep. Math. Phys., 5 (1), 1974, 121–130. [11] Meyer, K. R., Symmetries and integrals in mechanics, in Dynamical Systems, Peixoto, M. M. ed., Academic Press, New York, 1973, pp. 259–272. [12] Meyer, K. R., Periodic Solutions of the N –Body Problem, Lecture Notes in Math. 1719, Springer, Berlin, 1999. [13] Meyer, K. R. and Offin, D. C., Introduction to Hamiltonian Dynamical Systems and the N -Body Problem, Third Edition, Applied Mathematical Sciences 90, Springer, New York, 2017. [14] Meyer, K. R., Palacián, J. F. and Yanguas, P., Geometric averaging of Hamiltonian systems: periodic solutions, stability, and KAM tori, SIAM J. Appl. Dyn. Syst., 10 (3), 2011, 817–856. [15] Meyer, K. R., Palacián, J. F. and Yanguas, P., Stability of a Hamiltonian system in a limiting case, Regul. Chaotic Dyn., 17 (1), 2012, 24–35. [16] Meyer, K. R., Palacián, J. F. and Yanguas, P., Singular reduction of resonant Hamiltonians, Nonlinearity, 31 (6), 2018, 2854–2894. 17 [17] Moser, J., Regularization of Kepler’s problem and the averaging method on a manifold, Comm. Pure Appl. Math., XXIII, 1970, 609–636. [18] Reeb, G., Sur certaines propriétés topologiques des trajectoires des systèmes dynamiques, Acad. Roy. Belgique. Cl. Sci. Mém. Coll. in 8◦ , 27, 1952. [19] Robutel, P., Niederman, L. and Pousse, A., Rigorous treatment of the averaging process for co-orbital motions in the planetary problem, Comp. Appl. Math., 35 (3), 2016, 675–699. [20] Spirig, F. and Waldvogel, J., The three-body problem with two small masses: a singular- perturbation approach to the problem of Saturn’s coorbiting satellites, in Stability of the Solar System and Its Minor Natural and Artificial Bodies, Szebehely, V. G. ed., Reidel, Boston, MA, 1985, pp. 253–263. [21] Yanguas, P., Palacián, J. F., Meyer, K. R. and Dumas, H. S., Periodic solutions in Hamiltonian systems, averaging, and the Lunar problem, SIAM J. Appl. Dyn. Syst., 7 (2), 2008, 311–340. 18