This is a preprint of: “Analytical and numerical results on families of n-ejection-collision orbits in the RTBP”, Mercè Ollé, Òscar Rodrı́guez, Jaume Soler, Commun. Nonlinear Sci. Numer. Simul., vol. 90, 105294, 27 pp., 2020. DOI: [10.1016/j.cnsns.2020.105294] Analytical and numerical results on families of n-Ejection-collision orbits in the RTBP M. Ollé, O. Rodríguez, J. Soler Abstract. In the planar RTBP with mass ratio µ we regularise the singularity at one of the primaries by means of Levi-Civita’s transformation in a rotating frame. We solve the variational equations in a neighbourhood of the ejection/collision orbits, giving analytic expressions for the first terms in µ of the convergent expansion for orbits with eccentricity e ' 1. For high enough values of the Jacobi constant C we give√analytic expressions for the coefficients of the above expansion in powers of the small parameter 1/ C and we prove the existence of four families of the so called n-ejection-collision (EC) orbits, that are orbits which eject from the primary and reach n relative maxima in the distance with the primary before finally colliding with it. Moreover, massive numerical explorations extending the analytical result for any value of the mass ratio and bigger ranges of C are also shown and discussed. 1 Introduction We recall that the planar restricted three-body problem (PRTBP) is the problem of the motion of a massless particle in the gravitational field created by two bodies, the primaries, of mass µ (the mass ratio) and 1 − µ, which move around circular orbits around their centre of mass. Although this problem has been studied by many different authors, the so called ejection-collision orbits remain to be well understood. In this paper we state and prove a theorem on the existence of only four families of n-ejection- collision orbits, for any n ≥ 1, µ small enough and sufficiently large values of the Jacobi constant C. The argument that guarantees the existence of exactly four families relies on the application of a perturbative approach and the implicit function theorem as well. These families will be labelled by αn , βn , δn and γn . The case n = 1 was already proved in [5] using McGehee’s regularisation, which essentially consists of blowing up the singularity to an invariant manifold with two submanifolds of unstable equilibrium points which have stable and unstable manifolds. Ejection-collision orbits are then seen as heteroclinic connections between equilibrium points representing respectively the ejection and the collision. For µ = 0, every ejection orbit is a collision orbit after reaching its apocentre. For µ 6= 0 but small it was seen in [5] that only four of them survive. This approach, however, seems difficult to generalise to a class of orbits which eject from the primary, have n ≥ 1 apocentres with n−1 close approaches to the primary with no collision and finally collide with it (the n-EC orbits). The difficulty is due to the fact that it is not easy to follow, even numerically, an orbit which gets several times very close to a hyperbolic equilibrium point as the interval of time becomes unbounded. The way out has been to use Levi-Civita’s regularisation, which transforms the singularity at the primary into a regular point of a differential system, halfway between reaching the singularity unbounded velocity and completely stopping the motion there. The price to be paid is a double covering of the phase space, actually not a great nuisance, and a far greater complexity of the equations of motion, especially in a rotating frame. As far as we know, no analytical results are available regarding this kind of solutions. For details of the Levi-Civita regularization in the RTBP, see for example [4], [24] 1 M. Ollé, Ó. Rodríguez, J. Soler and very recently [20] (the author uses singular collision orbits associated with the second primary to find trajectories reaching the vicinity of the secondary, to low energies). So a main contribution of this paper is the generalization of the existence of four families of n-EC orbits for any n ≥ 1. The proof is based on considering Levi-Civita regularization and a perturbation approach. We remark that the computation of closed expressions for the solution of the variational equations for µ 6= 0 in a neighbourhood of the ejection-collision orbits is quite demanding and leads of course to integrals not expressible in closed form with elementary functions. However, the integrals √ depend on the Jacobi constant C and can be expanded in powers of the small parameter 1/ C, now with coefficients given by closed form expressions. Sections 5 and 6 are devoted to the actual computation of the expressions for the coefficients. If general values of µ and C are to be considered, the problem can on only be tackled numerically. This is done in section 7. Through numerical simulations, we have extended the previous analytical results, valid only for sufficiently small µ and C big enough, to the case all µ ∈ (0, 1) and large ranges of the Jacobi constant C (larger Hill’s regions). Given n ≥ 1, we show the existence of the four families αn , βn , δn and γn of n-EC orbits, for any value of the mass parameter µ ∈ (0, 1) and big ranges of values of C (we remove the restriction of C being necessarily large enough). We remark that besides illustrating the existence of four families of n-EC orbits to the big primary (that is µ ∈ (0, 0.5]), we also show the existence of four families of n-EC orbits to the small one (that is µ ∈ [0.5, 1)). Moreover, according to the theorem, given n and µ > 0 there is a value of Ĉ for which there exist exactly four families of n-EC orbits when C > Ĉ. We have computed this Ĉ value that provides a frontier between the region where there are exactly four families of EC orbits and the region where bifurcations of EC orbits take place and other families of EC orbits appear. Previous analytical results related to n-EC orbits in the planar circular RTBP are available only for the n = 1 case. We mention the paper by Llibre [13] where he proved the existence of at least two EC orbits for µ > 0 small enough and the Jacobi constant C big enough. The extension to the existence of four EC orbits for any µ ∈ (0, 0.5] and C big enough was done by Chenciner and Llibre [5]. Lacomba and Llibre [12] proved that both Hill problem and the RTBP have no C 1 -extensible regular integrals (taking into account EC orbits). Finally also some partial results related to the existence of 1-EC orbits of the spatial RTBP are given in [14] and of the planar elliptic RTBP (for µ and the eccentricity of the elliptical orbits of the primaries both positive and small enough) in [15] and [21]. We emphasize that the procedures used in the mentioned papers are mainly based on blow-up techniques (we refer to [7] and [16] for basic ideas) and the regularization of the equations of motion (see [23] and [24] for general theory in celestial mechanics). On the numerical side, we first mention Hénon’s paper [9], where the computation of EC orbits is done along the continuation of some families of symmetric periodic, non collision, orbits in the Copenhagen problem (that is µ = 0.5), in Hill’s problem (see [10]) and similarly in [2] the computation of 16 particular collision periodic orbits of the RTBP is done for various values of µ ∈ (0, 0.5]. As far as we know, the most complete papers about the numerical computation of n-EC orbits, for n = 1, ...., 25, the continuation of families and appearance of bifurcation orbits in the RTBP are [18] and [19]. We also mention the analysis of n-EC orbits in other contexts, for example see [17] in the atomic physics or in [1] in the collinear four-body problem. The paper is organized as follows: In Section 2, a short summary of the main properties of the RTBP are recalled, and in particular the equations of motion in the usual synodical (rotating) coordinates (with two singularities associated with collision with each primary) and in the rotating Levi-Civita ones (where the collision with the big primary has been regularized) are provided. In Section 3, the comparison between McGehee regularization and Levi-Civita one is shortly discussed (see more details in [19]) and the main geometrical ideas involved in the proof of the existence of n-EC orbits are presented. Section 4 contains the statement of the theorem and the core of the proof. Some detailed computations of the proof are done in Section 5 for the non perturbed 2-body problem (µ = 0) and in 2 Analytical and numerical results on families of n-Ejection-collision orbits in the RTBP Section 6 for the perturbed one (µ > 0 and small enough). We remark that our purpose has been to provide all the proofs (some of them maybe easy) of the Lemmas used to ease the readability of the paper and to make it self contained. Finally Section 7 is devoted to numerical simulations to extend the analytical results of the main theorem. We remark that the numerical integrations of the systems of ODE done along the paper use an own implemented Runge-Kutta (7)8 method with an adaptive step size control described in [8] (and a Taylor method implemented on a robust, fast and accurate software package by Jorba and Zou [11].) 2 The planar RTBP and the LC-transformation In this Section we shortly recall the circular, restricted three-body problem (RTBP) as well as some well known properties that will be used along the paper. In the RTBP we take two massive bodies P1 and P2 , called primaries, that describe circular orbits around their common center of mass, located at the origin. We consider a particle P with infinitesimal mass that moves on the same plane as the primaries under their gravitational forces. It will be useful to consider a coordinate system (called rotating or synodical system) of coordinates that rotates with the primaries, in such a way, that for suitable units of length, mass and time, the primaries will have mass 1 − µ and µ, µ ∈ (0, 0.5], their positions will be (µ, 0) and (µ − 1, 0) respectively, and the period of their motion will be 2π. In such context, the equations of motion for the particle in the rotating system are given by ( ẍ − 2ẏ = Ωx (x, y) (1) ÿ + 2ẋ = Ωy (x, y), where ˙ = d/dt and 1 2 1−µ µ 1 (x + y 2 ) + p Ω(x, y) = +p + µ(1 − µ) 2 2 (x − µ) + y 2 2 (x − µ + 1) + y 2 2 (2) 1 2 2  1 − µ µ = (1 − µ)r1 + µr2 + + 2 r1 r2 p p with r1 = (x − µ)2 + y 2 and r2 = (x − µ + 1)2 + y 2 . So, the equations become singular when r1 or r2 → 0. The main properties of this sytem used later on are the following (see [24] for details): 1. There exists a first integral, defined by C = 2Ω(x, y) − ẋ2 − ẏ 2 (3) and known as Jacobi integral. 2. System (1) has the symmetry (t, x, y, ẋ, ẏ) → (−t, x, −y, −ẋ, ẏ). (4) A geometrical interpretation of it is that given an orbit in the configuration space (x, y), the symmetrical orbit with respect to the x axis will also exist. 3. The simplest solutions are 5 equilibrium points: the so called collinear ones Li , i = 1, 2, 3, located on the x and the triangular ones Li , i = 4, 5. On the plane (x, y), L1,2,3 are located on the x axis, with x(L2 ) < µ − 1 < x(L1 ) < µ < x(L3 ) and L4,5 forming an equilateral triangle with the primaries. CLi will stand for the value of C at Li , i = 1, ..., 5. 3 M. Ollé, Ó. Rodríguez, J. Soler 4. Depending on the value of the Jacobi constant C, the particle can move on specific regions of the plane (x, y), called Hill regions and defined by  R(C) = (x, y) ∈ R2 | 2Ω(x, y) ≥ C . (5) In order to deal with the singularity of the first primary (r1 = 0) we will consider the Levi-Civita regularization (see [24]). The well known transformation of coordinates and time is given by:   x = µ + u2 − v 2   y = 2uv (6)   dt = 4 u2 + v 2   ds and we recall that each point in (x, y) coordinates is mapped to two points in (u, v) coordinates (see Figure 1). The system (1) becomes:  00    u − 8(u2 + v 2 )v 0 = 4U (u2 + v 2 ) u     8µu 8µu(u2 + v 2 )(u2 + v 2 + 1)   = 4µu + 16µu3 + 12(u2 + v 2 2 ) u + − − 4Cu     r2 r23  v 00 + 8(u2 + v 2 )u0 = 4U (u2 + v 2 ) v (7)     2 2 2 2 8µv 8µv(u + v )(u + v − 1)   = 4µv − 16µv 3 + 12(u2 + v 2 )2 v + − − 4Cv     r2 r23   0 C =0 where 0 = d/ds and 1h 2 i 1−µ µ C U= (1 − µ) u2 + v 2 + µ (1 + u2 − v 2 )2 + 4u2 v 2 + 2 + − . 2 u + v2 r2 2 p with r2 = (1 + u2 − v 2 )2 + 4u2 v 2 . The system of ODE (7) is now regular everywhere except at the collision with the small primary (r2 = 0). We remark that the system of ODE (7) makes sense for each value of a Jacobi constant C fixed. So, in order to take an initial condition of this system, we will take (u(0), v(0), u0 (0), v 0 (0), C(0)). However, for all the coming computations we will actually consider system (7) removing the last equation in C, and we will consider the corresponding solution for a fixed C and initial condition simply (u(0), v(0), u0 (0), v 0 (0)). In this system of reference, the previous properties of the RTBP are written as: 1. Jacobi Integral:  u02 + v 02 = 8 u2 + v 2 U (8) which is regular at the collision with the big primary. In particular (see [24]), the velocity at the position of the first primary (u = 0, v = 0) is given by: u02 + v 02 = 8(1 − µ) (9) 2. As the Levi-Civita transformation duplicates the configuration space (see Figure 1) the equations of motion satisfy two symmetries, (10a) in consequence of the duplication of space and (10b) due to (4): (s, u, v, u0 , v 0 ) → (s, −u, −v, −u0 , −v 0 ) (10a) 0 0 0 0 (s, u, v, u , v ) → (−s, −u, v, u , −v ) (10b) 4 Analytical and numerical results on families of n-Ejection-collision orbits in the RTBP Figure 1: Levi-Civita transformation. Hill’s region for µ = 0.2 and C = 3.53. Left. Synodic (x, y) coordinates. Right. Levi-Civita ones (u, v). In grey the complementary of Hill’s region. 3. The equilibrium points are now duplicated and they are located on the plane (u, v). In particular, the collinear points now are located in the u axis and in the v axis. 4. Similarly, given a value of the Jacobi constant C, the Hill’s region in variables (u, v) now becomes  R(C) = (u, v) ∈ R2 | (u2 + v 2 )U ≥ 0 . (11) In particular we will consider values of the Jacobi constant C ≥ CL1 , the value of the Jacobi constant associated to the equilibrium point L1 . In this way, it will be enough to regularize only the position of P1 because the Hill’s region associated to these values of C avoids collisions with the second primary (see Figure 2) Figure 2: Hill regions associated to CL1 (Left) and to C < CL1 (Right) for µ = 0.2. 5 M. Ollé, Ó. Rodríguez, J. Soler 3 n ejection-collision orbits. McGehee vs Levi-Civita As stated in the Introduction, the main goal of this paper is to prove analytically the existence of four families of n-EC orbits, for µ > 0 small enough and C big enough. We recall that an n-ejection- collision orbit, simply noted by n-EC orbit, is an orbit that ejects from the big primary and reaches n times a relative maximum in the distance to the big primary before colliding with it. Along this Section we introduce the main geometrical ideas and definitions that we will need for the proof of this result. As a motivation, we start this Section considering McGehee regularization and shortly recalling some results. We will also compare it with Levi-Civita one, and we explain why we will use Levi-Civita coordinates throughout the paper. Figure 3: (Left) Collision Manifold obtained with the McGehee regularization for any µ > 0. (Right) Ejection orbits manifold (in red) and collision orbits manifold (in blue) for µ = 0 and C > 0, in this case both manifolds are the same. The McGehee regularization consists in a blow up of the collision. The variables used are polar coordinates (r, θ) with respect to the primary to be regularized as the origin and the new variables dt , U = r V = r1/2 dr dt (not to be confused with potential U ) and a change of time dt/dτ = r are 3/2 dθ 3/2 introduced. The system of ODE in these variables is well defined at the origin (see [18]) and has an invariant torus Λ defined by r = 0, called the collision manifold (see Figure 3). On Λ there exist two circles of equilibrium points S + and S − characterized as follows: • Each equilibrium point P ∈ S + has an associated 2-d unstable manifold W u (P ) and a 1-d stable one W s (P ). • Each equilibrium point Q ∈ S − has an associated 2-d stable manifold W s (Q) and a 1-d unstable one W u (Q). The dynamics of the stable manifold of the points of S + and the unstable manifold of the points of S − corresponds to the internal dynamics on Λ due to the blow up of the collision that is made. The dynamics of the other manifolds associated with these equilibrium points are those that will allow us 6 Analytical and numerical results on families of n-Ejection-collision orbits in the RTBP to find the EC orbits. Actually, we will consider the ejection (collision) orbits for a fixed value of the Jacobi constant C, so the manifolds W u (P ) and W s (Q) are 1-d, and each ejection (collision) orbit can be characterized by a point P ∈ S + (Q ∈ S − ) or an angle θ ∈ [0, 2π). Thus, fixed C we have: • The ejection orbits manifold is the set of orbits on the unstable manifold W u (P ), when varying P ∈ S + . In this way, each ejection orbit corresponds to an orbit such that r > 0 for all finite time τ and tends to an equilibrium point P ∈ S + as τ → −∞. • The collision orbits manifold is the set of orbits on the stable manifold W s (Q), when varying Q ∈ S − . In this way, each collision orbit corresponds to an orbit such that r > 0 for all finite time τ and tends to an equilibrium point Q ∈ S − as τ → ∞. Hence each EC orbit corresponds to an heteroclinic connection between a point P ∈ S + and a point Q ∈ S − . And if we want to compute the heteroclinic connections between S + and S − we can do it via the intersection of the associated manifolds (see [18] for details). It should be noted that when the mass parameter µ is zero (and C > 0) the ejection orbits manifold and the collision one coincide, i.e. all the orbits that eject from the primary collide with it and, therefore, all the ejection orbits are 1-EC orbits. We show both manifolds in variables (V, r, θ) in Figure 3 right, where each circle S + and S − collapses to a point. But when µ 6= 0 both manifolds are different and the existence and number of n-EC orbits is no longer evident. In [18] the numerical computation and continuation (and bifurcation) of families of n-EC orbits were described and, since the McGehee’s coordinates were used there, any n-EC orbit was regarded as an heteroclinic orbit. At this point it is worthwhile comparing the two regularizations (McGehee and Levi-Civita) when applied to the study of ejection/collision orbits as a motivation to use Levi-Civita one throughout the paper. Ejection/collision orbits in the Levi-Civita regularization are simply orbits that leave from/arrive at the origin, which is now a regular point, so it takes a finite range of time to describe an EC orbit. By contrast, it takes an infinite time to describe an EC orbit in McGehee coordinates, since they are asymptotic (heteroclinic) connections. From this point of view, although the system of ODE in Levi-Civita variables is more intricated, the numerical computations are really faster. Moreover, the initial conditions of an ejection orbit are on invariant manifolds of equilibrium points when using McGehee variables, and (for numerical simulations) we take actually an approximation of such initial conditions. Conversely, in Levi-Civita p variables, we simply take initial conditions u = v = 0 and velocity any vector with norm 8(1 − µ) for a given C. Thus, from now on, we will only use Levi-Civita coordinates and the following definition. Definition 3.a. Given a value of the mass parameter µ, we define the ejection (collision) orbits manifold of energy level C as the set of orbits that have initial conditions p p (0, 0, 2 2(1 − µ) cos (θ0 ) , 2 2(1 − µ) sin (θ0 )), θ0 ∈ [0, 2π) integrated forward (backward) in time. In order to obtain the n-EC orbits we will use an idea explained in [18] for their numerical computation. Given a value of µ and C we compute the EC orbits as the intersection of the ejection orbits manifold and the collision orbits one. To do so we define the Poincaré section: Σ : {g(u) = uu0 + vv 0 = 0, g 0 (u) < 0}. (12) where u := (u, v, u0 , v 0 ). It should be noted that for every ejection or collision orbit, the crossing with this section corresponds to a maximum in the distance to the first primary. In this way, we define Dk+ as the k-th intersection of the ejection orbits manifold with Σ and, similarly we define Dk− as the k-th 7 M. Ollé, Ó. Rodríguez, J. Soler intersection of the collision orbits manifold with Σ. For µ > 0 small enough and C big enough, Dk+ and Dk− will be closed curves parametrized by an angle θ0 ∈ [0, 2π). In order to calculate the n-EC orbits we will use the following lemma: Lemma 1. The number of intersections between Di+ and Dj− corresponds to twice the number of (i + j − 1)-EC orbits. Proof. If an orbit belongs to the ejection orbits manifold and the collision orbits manifold, it is an ejection-collision orbit. In addition, the number of relative maxima in the distance to the first primary of such EC orbit is equal to the number of times that the orbit will cross Σ and that is i (following the ejection orbits manifold) plus j − 1 (following the collision orbits manifold). Finally, since the Levi-Civita regularisation doubles the configuration space, we have that the cardinal |Di+ ∩ Dj− | is twice the number of (i + j − 1)-EC orbits. Before enunciating the main Theorem it is important to observe that we do not need to compute Dk− , it can be obtained directly from Dk+ via the symmetry (10a) and for this reason we only work from now on with Dk+ . So a crutial point is to find the intersection points between the closed curves Di+ and Dj− . Any such intersection point gives rise forward (backward) in time to collision (ejection). Actually, for C given, these are curves in R4 . However, since by definition they are obtained as curves in Σ (this means two restrictions on the variables), we can consider only their projection in the (u, v) configuration space. In particular we can express Dk+ in polar coordinates, i.e. we will write, abusing 2 notation, Dk+ = (Rk+ (θ), θ), as a curve in R2 where Rk+ = (u2 + v 2 )|sk , sk being the necessary time to reach the k-intersection with Σ. Due to the expression that Rk has, we will work with the curve 2 2 2 Dk+ := (Rk2 (θ), θ). Thus, the intersections of Di+ and Dj− will be the same as the intersections of Di+ and Dj− . 4 Main Theorem The existence, the number and the characteristics of the n-ejection-collision orbits for small enough values of the mass parameter and for sufficiently restricted Hill regions can be summarized in the following Theorem. Theorem 1. For C big enough and for all n ∈ N there exists a µ̂(C, n) such that for µ ≤ µ̂(C, n) there exist four n-EC orbits, which can be characterized by: • Two n-EC orbits symmetric with respect to the x axis. • Two n-EC orbits with symmetric trajectory one of the other with respect to the x axis. The respective families (when varying C) are labelled by αn , γn , βn and δn . Proof. The proof is based on a perturbative approach. First of all, we can rewrite the equations of 8 Analytical and numerical results on families of n-Ejection-collision orbits in the RTBP motion (7) as:   0  0   u u0  0   2 2 2 2  v  v 0   8u 8u(u + v )(u + v + 1)   0 =    3  u   8(u2 + v 2 )v 0 + 12(u2 + v 2 )2 u − 4Cu  + µ 4u + 16u + r2 − r 3   2  v0 −8(u2 + v 2 )u0 + 12(u2 + v 2 )2 v − 4Cv  8v 8v(u2 + v 2 )(u 2 + v 2 − 1)  3 4v − 16v + − 3 r2 r2 = F0 (u, v, u0 , v 0 ) + µF1 (u, v) (13) and we look for n-EC orbits expressed in a series expansion in µ: u = u0 + µu1 + O(µ2 ) where u = u(s), u0 = u0 (s) and u1 = u1 (s). Developing with respect to µ we have: 0   u0 = u0 + µu1 + O(µ2 ) = F0 u0 + µu1 + O(µ2 ) + µF1 u0 + O(µ) (14) = F0 (u0 ) + µDF0 (u0 )u1 + µF1 (u0 ) + O(µ2 ) Therefore, u0 and u1 must satisfy 0 u0 = F0 (u0 ) (15) 10 0 1 0 u = DF0 (u )u + F1 (u ) (16) that is, u0 is a solution of the 2-body problem in Levi-Civita rotating coordinates with µ = 0 and u1 (s) is obtained as Z s u1 (s) = X(s)u10 + X(s) X −1 (s)F1 (u0 (s)) ds (17) 0 ! ∂u0 ∂ ∂ ∂ ∂ ∂ where X(s) := (s) and = , , 0 , 0 , being u00 = u0 (0) the initial conditions. ∂u00 ∂u00 ∂u00 ∂v00 ∂u00 ∂v00 The explicit computation of u0 (s) and u1 (s) will be done in Sections 5 and 6. Regarding the initial conditions, an ejection orbit will satisfy (u(0), v(0)) = (0, 0), and due to (9) the velocity is a vector with norm 8(1 − µ) and an initial angle θ0 ∈ [0, 2π) such that p p u0 = u(0) = (0, 0, 2 2(1 − µ) cos θ0 , 2 2(1 − µ) sin θ0 ) = u0 (0)+µu1 (0)+O(µ2 ) = u00 +µu10 +O(µ2 ) (18) with √ √ u00 = (0, 0, 2 2 cos θ0 , 2 2 sin θ0 ) √ √ (19) u10 = (0, 0, − 2 cos θ0 , − 2 sin θ0 ). As explained at the end of Section 3, in order to prove the theorem, we must compute the intersection 2 2 points between the curves Di+ and Dj− in R2 with i + j = n + 1. The precise number of intersection 2 points corresponds to the precise number of n-EC orbits. Moreover, once Dk+ is obtained, for some 2 k, Dk− is obtained applying the symmetry (10a). So a main goal is to compute the curve Dk+ , that is, for every initial condition parametrized by θ0 we want to compute the time sk necessary for the solution u with initial condition u0 to arrive at the k-th maximum distance (or the square of the distance). To do so we consider sk (θ0 ) = s0k + µs1k (θ0 ) + O(µ2 ) (20) 9 M. Ollé, Ó. Rodríguez, J. Soler where s0k (independent of θ0 ) is the known time for the 2-body solution u0 to reach the k-th maximum distance (the k-th crossing with Σ). In particular: (2k − 1)π Lemma 2. s0k = √ 4 C Proof. The formula is obtained directly from solving the unperturbed system (15) and imposing a maximum in the distance (see Section 5 for details). Then, via the Implicit Function Theorem we have: 0 0 0 0 u0 u1 + u0 u1 + v 0 v 1 + v 0 v 1 s1k (θ0 ) = − (21) u 0 02 +v 0 02 + u0 u0 00 + v 0 v 0 00 s0k Therefore, the square of the k-th maximum distance is given by:   2  2  u2 + v 2 sk (θ0 ) = u0 + µu1 + O(µ2 ) + v 0 + µv 1 + O(µ2 ) s0k +µs1k (θ0 )+O(µ2 ) " #      02 02 0 1 0 1 0 00 0 00 = u +v s0k + 2µ u u +v v s0k + s1k (θ0 ) u u +v v + O(µ2 ) s0k (22) In consequence, we have the k-th maximum distance (or square distance) in terms of θ0 and then work in polar coordinates, i.e. we can express Dk+ in polar coordinates as (Rk (θ0 ), θk (θ0 )) where Rk (θ0 ) is given by: q Rk (θ0 ) = (u2 + v 2 ) s (θ0 ) k q (23) 2 2 = Rk0 (θ0 ) + 2Rk1 (θ0 )µ + O(µ2 ) where 2  2 2  Rk0 (θ0 ) = u0 + v 0 s0k (24a) 2   0 0  Rk1 (θ0 ) = u0 u1 + v 0 v 1 s0 + s1k (θ0 ) u0 u0 + v 0 v 0 (24b) k s0k and θk (θ0 ) = θk0 (θ0 ) + O(µ). Also, thanks to the Poincaré section Σ that we are considering, in particular at the time s0k , we have  0 0  u0 u0 + v 0 v 0 =0 s0k Therefore it is not necessary to calculate the time s1k and the expression (24b) becomes: 2  Rk1 (θ0 ) = u0 u1 + v 0 v 1 s0k . (25) 2 2 The expressions of Rk0 (θ0 ) and Rk1 (θ0 ) are obtained in the following lemmas and the proofs of these can be found in Sections 5 and 6. 10 Analytical and numerical results on families of n-Ejection-collision orbits in the RTBP 2 2 Lemma 3. Rk0 = C Proof. The formula is obtained from solving the unperturbed system (15) and evaluating this solution at the time s0k (see Section 5 for details).  12 1 3 2 (3 cos(4θ0 ) + 1) 8 cos(2θ0 ) 3 − 5 cos2 (2θ0 ) Lemma 4. Rk (θ0 ) = − + 2 + + C C C4 C5 2(2k − 1)π sin(4θ0 ) + + O(C −13/2 ) C 11/2 Proof. The formula is obtained from solving the perturbed system (16) (see Section 6 for details). In addition, the expression of θk (θ0 ) is given by: (2k − 1)π Lemma 5. θk (θ0 ) = θ0 − + O(µ) 2C 3/2 Proof. The formula is obtained from the intrinsic rotation of the synodic system (see Section 6 for details). 2 Hence, the expression of Dk+ (θ0 ) in polar coordinates is 2  2 2  Dk+ (θ0 ) = Rk0 + 2µRk1 (θ0 ) + O(µ2 ) , θk (θ0 ) (26) and using the symmetry (10) we have: ! 2 2 2 (2k − 1)π Dk− (θ0 ) = Rk0 + 2µRk1 (θ0 ) 2 + O(µ ) , θ0 + + O(µ) (27) 2C 3/2 Finally, in order to compute the number of n-EC orbits it only remains to calculate the number of 2 2 intersections of Di+ and Dj− with i + j = n + 1: Lemma 6. For sufficiently small values of the mass parameter and n such that i + j = n + 1 there 2 2 exists a Ĉ(µ, n) such that for C ≥ Ĉ(µ, n) the curves Di+ and Dj− intersect at the 8 points: πk θ= + O(C −1 , µ), k = 0, ..., 7 (28) 4 2 2 Proof. Parametrizing the curves Di+ and Dj− with the same angle variable we have: ! ! 2 (2i − 1)π Di+ = Ri2 θ+ + O(µ) , θ (29) 2C 3/2 ! ! 2 (2j − 1)π Dj− = Rj2 θ− + O(µ) θ (30) 2C 3/2 with 11 M. Ollé, Ó. Rodríguez, J. Soler ! ! (2i − 1)π 2 2 (2i − 1)π Ri2 θ+ 3/2 + O(µ) = Ri0 + 2µRi1 θ+ + O(µ2 ) 2C 2C 3/2 ! ! (2j − 1)π 2 2 (2j − 1)π Rj2 θ− + O(µ) = Rj0 + 2µRj1 θ− + O(µ2 ) 2C 3/2 2C 3/2 2 2 In this way we can subtract the radii of Di+ and Dj− ! ! (2i − 1)π (2j − 1)π Ri2 θ+ − Rj2 θ− 2C 3/2 2C 3/2 " ! !# (31) 2 (2i − 1)π 2 (2j − 1)π = 2µ Ri1 θ+ − Rj1 θ− 2 + O(µ ) 2C 3/2 2C 3/2 and we will look for the values of θ in which this difference is null. Dividing by µ we have: " ! !# 12 (2i − 1)π 12 (2j − 1)π 2 Ri θ+ − Rj θ− + O(µ) = 0 (32) 2C 3/2 2C 3/2 Applying the Implicit Function Theorem for sufficiently small µ, we just have to find zeros (values of θ as a function of C) of ! ! 12 (2i − 1)π 12 (2j − 1)π 0 = Ri θ+ − Rj θ− 2C 3/2 2C 3/2 " ! !# 6 2(2i − 1)π 2(2j − 1)π = 4 cos 4θ + − cos 4θ − C C 3/2 C 3/2 " ! !! 8 (2i − 1)π 2 (2i − 1)π (33) + 5 cos 2θ + 3 − 5 cos 2θ + C C 3/2 C 3/2 ! !! # (2j − 1)π (2j − 1)π − cos 2θ − 3 − 5 cos2 2θ − C 3/2 C 3/2 " ! !# 2π 2(2i − 1)π 2(2j − 1)π + 11/2 (2i − 1) sin 4θ + − (2j − 1) sin 4θ − + O(C −13/2 ) C C 3/2 C 3/2 Similarly, multiplying by C 4 , expanding in C −1/2 and applying again the Implicit Function Theorem for sufficiently small C −1/2 , we just have to find the zeros of 4π (5i + 7j − 6) − sin (4θ) + O(C −5/2 ) (34) C 3/2 Therefore, the zeros of (32) are given by: πk θ= + O(C −1 , µ), k = 0, ..., 7 (35) 4 12 Analytical and numerical results on families of n-Ejection-collision orbits in the RTBP Finally, it only remains to verify that the zeros are transversal in θ and this can be easily seen by deriving the expression with respect to θ. Hence, equation (34) has eight and only eight zeros in θ as a function of C −1/2 , and by the Implicit Function Theorem equation (32) has the same number of zeros as a function of C −1 and µ. We have 2 2 that Di+ and Dj− intersect exactly at eight points. Using Lemma 1 we can conclude that there exist four and only four n-EC orbits. In addition, using the two symmetries of the problem we can see that those that have angle intersection θk with k = 0, 2 (4, 6) correspond to symmetric n-EC orbits (the (x, y) projection with respect to the x axis), i.e. αn and γn and those that have angle intersection θk with k = 1, 3 (5, 7) correspond to orbits with symmetric trajectory one with respect to the other (also the (x, y) projection), i.e. βn and δn . Thus, it only remains to prove Lemmas 2, 3 and 5 (in Section 5) and Lemma 4 (in Section 6). In order to do it the first step is to solve the system (15). This will be done in the next Section. Remark 1. Since the proof of Theorem 1 is based on the application of the Implicit Function Theorem, the assertion states that given C big enough and any n, there exists a sufficient small value of µ > 0 for which there exist four n-EC orbits. As a theoretical result it is difficult to find such value of µ depending on C and n, and the value of µ must be small enough. But we are interested in obtaining n-EC orbits for any value of µ ∈ (0, 1). Therefore we apply a numerical strategy that, for any fixed value of µ ∈ (0, 1), enables us to show the existence of four n-EC orbits for big ranges of C. This will be our numerical approach in Section 7 devoted to numerical simulations. We will fix a value of µ and, given n, we will find a value Ĉ(µ, n) depending on µ and n such that for any C ≥ Ĉ(µ, n) there exist four n-EC orbits. 5 Non perturbed case (µ = 0) In order to find an explicit expression for the solution u0 (s) of the non perturbed system (µ = 0) (that is the 2-body problem) in rotating Levi-Civita (LC) coordinates, we will consider previously the problem in sidereal (non-rotating) coordinates, since it is simpler to solve and will be useful to obtain the first order of the perturbed solution u1 (s). In the next subsection we will prove Lemmas 2 and 3. 2-body problem in LC sidereal coordinates It is well known that the problem of two bodies (P1 and P ), where (X, Y ) is the position of P in sidereal coordinates is described by the system of ODE  2   d X X  2 =− 3 dt R (36)  2 d Y =−Y  dt2 R3 √ where R = X 2 + Y 2 and with a singularity at the origin (collision) (X, Y ) = (0, 0). The Levi-Civita transformation to new coordinates (û, v̂) and new time s defined by 13 M. Ollé, Ó. Rodríguez, J. Soler  2 2   X = û0 − v̂ 0   Y = 2û0 v̂ 0 (37)      dt = 4 û02 + v̂ 02  ds regularizes the ODE giving rise to the new system (see [24] for details)  00 û0 = −4K û0 00 (38) v̂ 0 = −4K v̂ 0 where 0 = d/ds, K = −2E and E is the energy integral. Taking into account the relationship between E and the Jacobi constant C and the angular momentum M , E = − C2 + M and the fact that 0 0 M = X Ẏ − Y Ẋ becomes M = (û0 v̂ 0 − v̂ 0 û0 )/2, the previous system becomes  000 0 0 û = −4(C − û0 v̂ 0 + v̂ 0 û0 )û0 00 0 0 (39) v̂ 0 = −4(C − û0 v̂ 0 + v̂ 0 û0 )v̂ 0 As it is well known the system (39) preserves the angular momentum, so we have: 0 0 Lemma 7. The quantity −û0 v̂ 0 + v̂ 0 û0 is constant along the trajectories 00 00 Proof. Simply by deriving and replacing the values of û0 and v̂ 0 we have: d 0 0 00 00 0 0 0 0 − û0 v̂ 0 + v̂ 0 û0 = −û0 v̂ 0 + v̂ 0 û0 = 4(C − û0 v̂ 0 + v̂ 0 û0 )û0 v̂ 0 − 4(C − û0 v̂ 0 + v̂ 0 û0 )û0 v̂ 0 = 0 ds Keeping that in mind, the solution of (39) with initial conditions û0 (0) = û00 is given by:  0   √ û00 √  0 0 û (s) = û0 cos(2 Ks) + √ sin(2 Ks) 2 K 0 (40)   √ v̂00 √   v̂ 0 (s) = v̂ 0 0 cos(2 Ks) + √ sin(2 Ks) 2 K 0 0 where K = C − û00 v̂00 + v̂00 û00 and has derivatives:  00 √ √ 0 √ û (s) = −2√K û00 sin(2√Ks) + û00 cos(2√ Ks) 0 0 (41) v̂ 0 (s) = −2 K v̂00 sin(2 Ks) + v̂00 cos(2 Ks) √ √ In particular, we are interested in ejection orbits, which have initial conditions (0, 0, 2 2 cos θ0 , 2 2 sin θ0 ), θ0 ∈ [0, 2π), i.e. √   2 cos θ0 √   û 0 (s) = √ sin(2 Cs) √ C (42)   2 sin θ0 √ v̂ (s) = √ 0 sin(2 Cs). C With this equation it is immediate to see the following two results, because the maximum distance does not vary with or without rotation: (2k − 1)π • s0k = √ . 4 C 14 Analytical and numerical results on families of n-Ejection-collision orbits in the RTBP  2  2 2 2 • Rk0 = û0 + v̂ 0 (s0k ) = . C which are the statements of Lemma 2 and Lemma 3 respectively. In order to prove Lemma 5, it is necessary to study the system in synodical coordinates with the purpose of calculating the rotation consequences (in particular the θk0 value) and because we will also need the 2-body solution in the synodical coordinates to calculate the first order solution u1 . 2-body problem in LC synodical coordinates In order to have the solution u0 (s) in Levi-Civita rotating frame, we simply apply the rotation transformation to the solution (40) in Levi-Civita sidereal coordinates. This is shown in the following lemma. Lemma 8. The solution of (15) is given by u0 (s) = û0 cos(−t/2) − v̂ 0 sin(−t/2) v 0 (s) = û0 sin(−t/2) + v̂ 0 cos(−t/2)  0   0  0 2 2 2 2 (43) u0 (s) = û0 + 2(û0 + v̂ 0 )v̂ 0 cos(−t/2) − v̂ 0 − 2(û0 + v̂ 0 )û0 sin(−t/2) 0  0 2 2   0 2 2  v 0 (s) = û0 + 2(û0 + v̂ 0 )v̂ 0 sin(−t/2) + v̂ 0 − 2(û0 + v̂ 0 )û0 cos(−t/2) 0 0 0 0 where û0 , v̂ 0 , û0 , vˆ0 and t stand for û0 (s), v̂ 0 (s), û0 (s), v̂ 0 (s) and t(s). Proof. We simply make use of the change of variables (using complex notation) • W = û + iv̂ 7−→ Z = X + iY by Z = W 2 , where (X, Y ) are the cartesian sidereal variables. • Z 7−→ z = x+iy by z = e−it Z, i.e. the rotation to transform from sidereal to synodic coordinates (x, y). • z 7−→ w = u + iv by z = w2 , being (u, v) the synodical Levi-Civita variables. t So the relation w = e− 2 i W for the 2-body problem can be written as u0 (s) = û0 cos(−t/2) − v̂ 0 sin(−t/2) (44) v 0 (s) = û0 sin(−t/2) + v̂ 0 cos(−t/2) dt On the other hand, u0 (s) and v 0 (s) are obtained using = 4(û2 + v̂ 2 ). ds The relation between the sidereal initial conditions and the synodical ones can be obtained from the previous lemma  0  0  u0 = û00   û0 = u0       v00 = v̂00   v̂00 = v0 0 0 2 2 0 0 2 2 (45)    u00 = û00 + 2(û00 + v̂00 )v̂00    û00 = u00 − 2(u00 + v00 )v00   00   0 2 2 v0 = v̂00 − 2(û00 + v̂00 )û00  v̂ 00 = v 0 0 + 2(u0 2 + v 0 2 )u0 0 0 0 0 0 15 M. Ollé, Ó. Rodríguez, J. Soler Remarks: 1. In the particular case of √ ejection √ orbits, the initial conditions both in sidereal and synodical coordinates are u00 = (0, 0, 2 2 cos θ0 , 2 2 sin θ0 ), θ0 ∈ [0, 2π). 2. We also notice that if we express the position of the particle with Levi-Civita sidereal variables, using polar complex notation, that is, W = Reiθsid , where θsid refers to the polar angle, then we have the corresponding two points in Levi-Civita synodical variables; using polar complex notation we obtain w+ = Rei(θsid − 2 ) t (46a) w− = Rei( θsid +π− 2t ) (46b) but we can omit expression (46b) due to the symmetry (10a). Now let us obtain an explicit expression for t(s). Lemma 9.  2 2   2 2  √ 0 û0 û0 + v̂ 0 v̂ 0 0 √ √ (a) u0 + v 0 (s) = û00 + v̂00 cos2 (2 Ks) + 0 0√ 0 0 cos(2 Ks) sin(2 Ks) (47) K 002 002 √ û + v̂0 + 0 sin2 (2 Ks) 4K 02 02 ! 02 02 ! √ 2 2 û0 + v̂00 2 2 û0 + v̂00 sin(4 Ks) (b) t(s) = 2(û00 + v̂00 ) + 0 s+ 4(û00 + v̂00 ) − 0 √ (48) 2K K 8 K 0 0 û00 û00 + v̂00 v̂00 √ + sin2 (2 Ks) K Proof. (a) It follows directly from (40). (b) It follows directly from the integration Z 2 2 t(s) = 4(û0 + û0 ) ds (49) In particular we have for the ejection orbits √ 4 sin(4 Cs) t(s) = s − . (50) C C 3/2 With this result and using (46), we can compute easily the term θk0 , i.e. the accumulated rotation during the span of time of length s0k t(s0k ) (2k − 1)π θk (θ0 ) = θ0 − + O = θ0 − +O (51) 2 2C 3/2 which is the statement of Lemma 5. It should be noted that to this rotation we should add an angle of π for even orbits, or equivalently add (k − 1)π to θk (θ0 ) (see (42)), but we can omit this increment thanks to the symmetry (10a) of the problem. 16 Analytical and numerical results on families of n-Ejection-collision orbits in the RTBP 6 Perturbed Problem 2 It only remains to prove Lemma 4. To do so, we need to compute u1 (s) and Rk+ (θ0 ). As previously mentioned, the solution of (16) is given by (17), therefore, the first step is compute X in a cunning way. To do so we will use the relationship between sidereal and synodic coordinates, so from (43) we have: ∂u0 (s) = V1 cos(−t/2) − V2 sin(−t/2) ∂u00 ∂v 0 (s) = V1 sin(−t/2) + V2 cos(−t/2) ∂u00 0 (52) ∂u0 (s) = V3 cos(−t/2) − V4 sin(−t/2) ∂u00 0 ∂v 0 (s) = V3 sin(−t/2) + V4 cos(−t/2) ∂u00 where:   ∂ û0 v̂ 0 ∂t V1 ≡ V1 (s) = + (s) ∂u00 2 ∂u00  0  ∂v̂ û0 dt V2 ≡ V2 (s) = − (s) ∂u00 2 ∂u00 " # ∂ û0 0  2  ∂v̂ 0 0 0 2 2 v̂ 0 − 2(û0 + v̂ 0 )û0 ∂t (53) 0 02 0 0 ∂ û V3 ≡ V3 (s) = + 2 û + 3v̂ + 4û v̂ + (s) ∂u00 ∂u00 ∂u00 2 ∂u00 " # ∂v̂ 0 0  2  0 0 0 2 2 û0 + 2(û0 + v̂ 0 )v̂ 0 ∂t 0 02 ∂ û 0 0 ∂v̂ V4 ≡ V4 (s) = − 2 3û + v̂ − 4û v̂ − (s) ∂u00 ∂u00 ∂u00 2 ∂u00 i.e.:    cos(−t/2) − sin(−t/2) 0 0 V1  sin(−t/2) cos(−t/2) 0 0  V2  X = RV =     (54) 0 0 cos(−t/2) − sin(−t/2) V3  0 0 sin(−t/2) cos(−t/2) V4 where V is the matrix that has rows given by Vi , i = 1, ..., 4. Thus, we only need to compute previously: ∂ û0 ∂ û0 ∂ û00 (s) = (s) ∂u00 ∂ û00 ∂u00 (55) ∂t ∂t ∂ û00 0 (s) = 0 (s) 0 ∂u0 ∂ û0 ∂u0 ∂ û00 The computation of is immediate from (45) and in the case of ejection orbits it is simply the ∂u00 identity√ matrix Id√since we just need to substitute in this matrix the values of the initial conditions (0, 0, 2 2 cos θ0 , 2 2 sin θ0 ) for an ejection orbit, 17 M. Ollé, Ó. Rodríguez, J. Soler   1 0 0 0 ∂ û00 √ √  0 1 0 0 (0, 0, 2 2 cos θ0 , 2 2 sin θ0 ) =   −4u00 v00 2 2  = Id. ∂u00 −2u00 − 6v00 1 0 2 2 6u00 + 2v00 4u00 v00 0 1 √ √ (0,0,2 2 cos θ0 ,2 2 sin θ0 ) (56) ∂ û0 ∂t To facilitate the reading of this section, the values of , and other values necessary for the ∂ û00 ∂ û00 computation of V can be found in the Appendix A. Similarly, the complete expression of the matrix V and V −1 can be found in it. Once the matrix V −1 has been calculated, we can compute the product X −1 F1 = V −1 R−1 F1 but the resulting expression to be integrated √ is quite unwieldy. The first step to simplify the expression is to introduce the variable τ = 2 Cs (see Appendix B). In this way we eliminate the terms related with C from most trigonometric expressions, except those that have a direct relationship with the rotation (i.e. with t(s)) √ sin(2 Cs) = sin(τ ) √ !   2 sin(4 Cs) 1 −2τ + sin(2τ ) sin(−t/2) = sin − s + = sin C 2C 3/2 C 3/2 2 and similarly with the cosine terms. Despite this change, the resulting expression cannot be integrated in finite form√ with elementary functions. To deal with this last expression we consider a series expansion in 1/ C. This expansion can be done since 1 √ • √ V −1 has a closed expression: constant terms and terms of higher order in 1/ C. 2 C • R−1 is simply a rotation. p • F1 depends only on u and v, and the denominator r1 is of the form 1 + O(C −1 ), therefore it can be expanded. The computation details are given in Appendix B. Thus, the integral in (17) becomes Z s Z τ      1 τ τ X −1 (s)F1 (u0 (s))ds = √ X −1 √ F1 u0 √ dτ. (57) 0 0 2 C 2 C 2 C The expressions  obtained after these considerations (simplifying the notation) are as follows, where X −1 F1 i is the i-th component of X −1 F1 √ √  1 −1  3 2 cos θ0 sin2 τ 24 2 cos θ0 2 cos4 θ0 − 1 sin6 τ √ X F1 1 = − − + O(C −9/2 ) 2 C C 3/2 C 7/2 √ √  1 −1  3 2 sin θ0 sin2 τ 24 2 sin θ0 2 sin4 θ0 − 1 sin6 τ √ X F1 2 = − − + O(C −9/2 ) 2 C C 3/2 C 7/2 18 Analytical and numerical results on families of n-Ejection-collision orbits in the RTBP √ √ 1  6 2 cos θ0 cos τ sin τ 12 2 sin θ0 sin τ (sin τ − τ cos τ ) √ X −1 F1 3 = + 2 C C C 5/2 √  48 2 cos θ0 2 cos θ0 − 1 cos τ sin5 τ 4 + √ C3 7 128 2 sin τ cos τ cos θ0 (5 cos6 θ0 − 6 cos2 θ0 + 2) − √ C4 384 2 h 4 3 + sin τ cos(τ /2) sin θ 0 32 cos (τ /2) C 9/2 + 52 cos3 (τ /2) sin2 θ0 + 30 cos3 (τ /2) sin4 θ0 + (28 − 64 sin2 θ0 + 32 sin4 θ0 ) cos5 (τ /2)(cos2 (τ /2) − 2)  i + −2 + 6 sin2 θ0 + sin4 θ0 (2 cos (τ /2) − τ sin (τ /2) cos τ ) √ 320 2 sin9 τ cos τ cos θ0 (14 cos8 θ0 − 30 cos4 θ0 + 20 cos2 θ0 − 3) + C5 + O(C −11/2 ) √ √ (58) 1 −1  6 2 sin θ0 cos τ sin τ 12 2 cos θ0 sin τ (sin τ − τ cos τ ) √ X F1 4 = − 2 C C C 5/2 √  48 2 sin θ0 2 sin4 θ0 − 1 cos τ sin5 τ + √ C3 128 2 sin τ cos τ sin θ0 (5 sin6 θ0 − 6 sin2 θ0 + 2) 7 − √ C4 384 2 h 4 − sin τ cos(τ /2) cos θ 0 32 cos3 (τ /2) C 9/2 + 52 cos3 (τ /2) cos2 θ0 + 30 cos3 (τ /2) cos4 θ0 + (28 − 64 cos2 θ0 + 32 cos4 θ0 ) cos5 (τ /2)(cos2 (τ /2) − 2)  i + −2 + 6 cos2 θ0 + cos4 θ0 (2 cos (τ /2) − τ sin (τ /2) cos τ ) √ 320 2 sin9 τ cos τ sin θ0 (14 sin8 θ0 − 30 sin4 θ0 + 20 sin2 θ0 − 3) + C5 −11/2 + O(C ) Integrating these results with respect to s on the interval [0, s0k ] (see Lemma 2) or [0, 2k−1 2 π] with respect to τ (see Appendix B) and using (62) we obtain the following expression for u1 (s0k ) √ √ 1 (−1)k 2 cos θ0 3(−1)k 2 cos θ0 u (s0k ) = √ − 2 C 2C 3/2 √ √ 3(−1)k (2k − 1) 2π sin θ0 3(−1)k (2k − 1) 2π sin θ0 + − √ 4C 2 2C 3 k 2 2 4 (−1) 2 cos θ0 (5π (2k − 1) + 128 cos θ0 − 64) − √ 16C 7/2 (−1) 2 cos θ0 (27(2k − 1)2 π 2 + 1920 cos6 θ0 − 2304 cos2 θ0 + 768) k + + O(C −6 ) 48C 9/2 19 M. Ollé, Ó. Rodríguez, J. Soler √ √ 1 (−1)k 2 sin θ0 3(−1)k 2 sin θ0 v (s0k ) = √ − 2 C 2C 3/2 k √ √ 3(−1) (2k − 1) 2π cos θ0 3(−1)k (2k − 1) 2π cos θ0 − + √ 4C 2 2C 3 (59) k 2 2 4 (−1) 2 sin θ0 (5π (2k − 1) + 128 sin θ0 − 64) − √ 16C 7/2 (−1) 2 sin θ0 (27(2k − 1)2 π 2 + 1920 sin6 θ0 − 2304 sin2 θ0 + 768) k − + O(C −6 ) 48C 9/2 2 Finally, from (25) we obtain the following expression for Rk1 (θ0 )  2 1 3 2(3 cos(4θ0 ) + 1) 8 cos(2θ0 ) 3 − 5 cos2 (2θ0 ) Rk1 (θ0 ) =− + 2 + + C C C4 C5 (60) 2(2k − 1)π sin(4θ0 ) + + O(C −13/2 ), C 11/2 which finishes the proof of Lemma 4. 7 Numerical results Once we have proved analytically, by a perturbation approach, the existence of exactly four families of n-EC orbits, for n given, µ small enough and for sufficiently large values of the Jacobi constant C, that is, C ≥ Ĉ (Theorem 1), the present Section is devoted to the extension of this result from a numerical point of view. More precisely, our purpose is twofold: (i) We have extended the existence of four families of n-EC orbits in both µ and C. On the one hand, we show that there exist exactly four families not only to the big primary, that is, for any µ ∈ (0, 0.5], but also the existence of exactly four families of n-EC orbits to the small one, i. e., for any µ ∈ [0.5, 1). Of course, from a perturbative point of view, µ must be small enough, but for numerical simulations, any value of µ can be considered. On the other hand, we have taken less restrictive values of C (not necessarily large enough), that is, bigger Hill’s regions. (ii) We have computed the Ĉ value according to the Remark of Theorem 1 that provides a frontier between the region where there are exactly four families of EC orbits and the region where bifurcations of EC orbits take place and other families of EC orbits appear. Figure 4: µ = 0.1, n = 2, C = 8. (Left) Curves D2+ (in red) and D1− (in blue) in the plane (θ, R). (Right) The four 2-EC orbits plotted in the synodic coordinates (x, y). 20 Analytical and numerical results on families of n-Ejection-collision orbits in the RTBP Concerning item (i), and recalling that a clue point for the existence of four n-EC orbits in the proof of Theorem 1 is the intersection of curves Di+ and Dj− (with i + j = n + 1) for µ ∈ (0, 1), n and C given, we have done the following computations. First, we take a value of µ (not small enough) and n fixed. We fix a big value of C and we have computed curves Di+ and Dj− such that i + j = n + 1, integrating the unstable ejection and stable collision manifolds, using the RTBP equations in Levi-Civita variables (that is, system (7)), up to the Poincaré section Σ. We have checked that both curves intersect in 8 transversal points (in rotating Levi-Civita coordinates), which correspond to four n-EC orbits (in cartesian rotating (x, y) ones). In Figure 4 left, curves D2+ and D1− (the projection in polar coordinates (R, θ)) are plotted for µ = 0.1, n = 2 and C = 8. So, four 2-EC orbits are obtained and plotted in Figure 4 right, in the synodical coordinates (x, y). Any such trajectory ejects from the big primary, has a maximum passage to it, a minimum one (with non zero distance), another maximum passage and finally collides with it. In Table 1 we provide the initial angle θ0 for √ each 2-EC orbit. √ So the associated initial conditions (in Levi-Civita synodical coordinates) are (0, 0, 2 1.8 cos θ0 , 2 1.8 sin θ0 ). α2 β2 γ2 δ2 θ0 1.7037895180958 2.44063260106480 0.1327163661903 0.9650603806965 Table 1: µ = 0.1, C = 8. Values of θ0 for 2-EC orbits. Figure 5: Curves D2+ and D1− for µ = 0.1 varying C. Second, we want to show the evolution of curves Di+ and Dj− , for µ (not small enough) and n given, when varying C, for a range of C (C big). We provide in Figure 5, the 3D plot corresponding to the (now) surfaces D2+ and D1− in the variables θ, R and C (for µ = 0.1). Therefore, for any C in a (small) given range, there exist the 8 intersecting points (in Levi-Civita variables) corresponding to the 4 n-EC orbits in synodical cartesian coordinates (x, y). So we can conclude, from our numerical simulations, that the four families of n-EC orbits persist taking bigger values of µ > 0 (not necessarily small enough) and a range of values of C (not necessarily big enough). Third, we have extended the previous results for any value of the mass parameter µ ∈ (0, 1) and big ranges of values of C. So, fixed any µ ∈ (0, 1) and for any n given, we have done the numerical 21 M. Ollé, Ó. Rodríguez, J. Soler continuation of the four families (labelled γn , δn , αn and βn ) of n-EC orbits in a big range of C. We plot in Figure 6 top left such families for µ = 0.1 and n = 1, 2, 5 (blue, green and purple colours respectively) in the plane (θ0 , C), where we recall that θ0 is the angle, in Levi-Civita polar coordinates, for the velocity in the initial condition of an ejection or collision orbit (see Definition 3). In particular for C = 5.5 we plot the corresponding n-EC orbits in Figure 6 top right in (x, y) variables. We emphasize that taking values of µ in the whole interval (0, 1), we are covering also the existence of four families of n-EC orbits with the small primary. As a particular example, we plot in Figure 6 middle left the four families of n-EC orbits, for n = 1, 2, 5 and µ = 0.9 in the (θ0 , C) plane. In Figure 6 middle right we plot the corresponding EC orbits for C = 5.5 and in Figure 6 bottom, we show the 3-EC orbits corresponding to each primary for µ = 0.4 and the same value of C = 4.5 (together with the Hill’s regions in (x, y) coordinates). We also remark that not only any value of µ can be considered but, of course, any value of n. For example we plot four n-EC orbits in Figure 7 top left, for n = 15, and top right for n = 25. It is worth mentioning two specific effects. First, given an n-EC orbit, let us define the interval [θe , θc ] that covers the range angle between the velocity vector at ejection θe and the one at collision θc . As far as n increases, since the required time to go from ejection to collision also does increase, for the same value of C, the range [θe , θc ] will also be larger for bigger n. This effect is clearly seen comparing Figure 6 for n = 5 with Figure 7 for n = 15 and n = 25. A second effect is on the shape of the petals of an n-EC orbit. It is clear that as far as µ increases and, more particularly, C decreases, the Hill’s regions are larger and the influence of both the small primary and the invariant objects that exist in such regions (and that did not exist for larger values of C, i.e. more restrictive Hill’s regions) is more visible. This effect can be observed comparing Figure 6 top and middle, where the petals look the same shape, versus the bottom Figure where we can appreciate a subtle difference in the shape of the three petals (the one in the middle is either thinner or fatter). Also in the 3D plot of Figure 7 bottom, we show the evolution of the 1-EC orbits (on the plane (x, y)) when varying C (on the vertical axis) for µ = 0.1; we can appreciate how the size of the petals gets bigger for decreasing C. This effect is shown for n = 1 for clearness but of course takes place for any value of n. Finally, we show in Figure 8 a 3D plot with the continuation of the four families of n-EC orbits, for n = 1, 2, 5, in the plane (θ0 , C) (θ0 ∈ [0, π] and C ∈ [5, 20]) when varying µ ∈ (0, 1) (in the vertical axis). From this plot it is worth making some remarks. When C is big enough, we see that the values of θ0 tend to πk/4, for k = 0, 1, 2, 3 for any value of µ. This is true as proved from Theorem 1, for small enough values of µ. A proof for any value of µ is left for a next future. At this point we mention references [6] and [12] for the proof of 1-EC orbits in the Hill’s problem, as a model for the motion close to the small primary. Another important comment is for small values of µ and when C decreases. Let us focus on the curves in the (θ0 , C) plane, for θ0 ∈ [0, π] and C ∈ [5, 20], and µ small fixed. When n = 1, such four curves look like straight lines; however as n increases, such lines become more curved when C decreases (that is when the Hill’s regions get bigger). If even smaller values of C are considered, the merging of families takes place and bifurcations of new families appear. These phenomena also happen for bigger values of µ. See [18] for details. Focussing on item (ii) and according to Theorem 1, for C big enough and n given, there exist values of µ small enough such that the existence of four n-EC orbits is proved. However, we want to consider any value of µ ∈ (0, 1). So, as explained in Remark 1, we will take a fixed value of µ ∈ (0, 1) and, given n, we will compute the value Ĉ(µ, n) such that for C ≥ Ĉ, there appear exactly four n-EC orbits. So for C < Ĉ we expect a bigger number of EC orbits, that is, bifurcations of EC orbits turn up. The numerical method to compute Ĉ, for a fixed µ ∈ (0, 1) and given n, is the following: we take an interval I = [Ca , Cb ] of values of C, and for each C ∈ I, (starting at Cb ) we vary θ0 ∈ [0, π) (that defines the initial conditions in synodical Levi-Civita variables) and find the four specific values of θ0 corresponding to four n-EC orbits. We say that we have four n-EC orbits for that C. So decreasing C we obtain four families of n-EC orbits. See Figure 6 top and middle left. However as we decrease C, we find a value of C ∈ I such that more than four n-EC orbits are found. This means that new families have bifurcated. Next we refine the value of C such that it is the frontier before appearing 22 Analytical and numerical results on families of n-Ejection-collision orbits in the RTBP Figure 6: Families γn , δn , αn and βn of n-EC orbits for n = 1, 2, 5 (blue, green and purple colours respectively) in the plane (θ0 , C): EC orbits to the big primary (top left) for µ = 0.1, and to the small one (middle left) for µ = 0.9. The corresponding n-EC orbit in (x, y) variables for C = 8 (top and middle right). (Bottom) 3-EC orbits to each primary for µ = 0.4 and the same value of C = 4.5 (together with the Hill’s regions). 23 M. Ollé, Ó. Rodríguez, J. Soler Figure 7: (Top) n-EC orbits to the big primary for µ = 0.1 in (x, y) coordinates. (Left) n = 15 and C = 15. (Right) n = 25 and C = 20. (Bottom) Evolution of the 1-EC orbits in (x, y) coordinates, varying C (the vertical axis). new families of n-EC orbits. That is precisely the specific values of Ĉ. In Figure 9 we show this kind of computation for µ = 0.1. We take the interval I = [CL1 , Cb ], Cb big enough. The reason to take Ca = CL1 is to guarantee that there is no possible path between P1 and P2 , due to the topology of the Hill region, so, roughly speaking, the dynamics in each bounded region around each primary is the simplest possible. In the figure, the x axis corresponds to θ0 ∈ [0, π) and the y one to C. The blue and red curves correspond to the families n = 2 and n = 3 EC orbits, respectively. We see on the plot the specific values Ĉ(0.1, 2) = 3.72442505 and Ĉ(0.1, 3) = 3.80644009. Clearly, for C ≥ Ĉ there exist four families of n-EC orbits whereas new families of n-EC orbits appear for C < Ĉ. For the fixed value of µ = 0.1, we have computed the Ĉ(0.1, n) value for each value of n varying from 24 Analytical and numerical results on families of n-Ejection-collision orbits in the RTBP Figure 8: Continuation of families γn , δn , αn and βn of n-EC orbits for n = 1, 2, 5 (blue, green and purple colours respectively) in θ0 and C when varying µ ∈ (0, 1) (on the vertical axis). 2 to 20. In Figure 10 we can see the evolution of Ĉ(0.1, n) (on the vertical axis) increasing n from 2 to 20 (on the horizontal one). We observe that the Ĉ(0.1, 1) is not on the plot because bifurcations appear when C < CL1 in this case. As expected, when n increases, a bigger value of Ĉ is obtained, that is, bifurcations of EC orbits appear earlier, in the sense that Ĉ is bigger. On the contrary, for smaller values of n, Ĉ is smaller; so we obtain larger ranges for C, meaning that for C ∈ (Ĉ, ∞) there exist four and only four families of n-EC orbits. We refer the interested reader to [18] for the details related to the bifurcation of families of EC orbits. 25 M. Ollé, Ó. Rodríguez, J. Soler Figure 9: µ = 0.1, families of n-EC orbits in this range of values of C (vertical axis), blue curves for n = 2 and red for n = 3. θ0 ∈ [0, π) is on the x axis. The specific values Ĉ(0.1, 2) and Ĉ(0.1, 3) as frontier of bifurcation families is also shown. Figure 10: For µ = 0.1, value of the frontier value Ĉ (on the y axis) versus n (on the x axis). 26 Analytical and numerical results on families of n-Ejection-collision orbits in the RTBP 8 Acknowledgments M. Ollé and O. Rodrıguez were supported by the Spanish MINECO/FEDER grant PGC2018-100928- B-100 and the Catalan grant 2017SGR-1049. J. Soler was supported by MINECO/FEDER grant number MTM2016-77278-P. References [1] Alvarez-Ramírez, M., Barrabés, E., Medina, M. and Ollé, M. Ejection-collision orbits in the symmetric collinear four-body problem. Commun Nonlinear Sci Numer Simulat 2019;71:82–100. [2] Bozis, G. Sets of collision periodic orbits in the Restricted problem. In: Periodic orbits, stability and resonances, G.E.O. Giacaglia (eds). Holland: D. Reidel Pub. Co.; 1970, p. 176–191. [3] Brunello, A. F., Uzer, T., Farrelly, D. Hydrogen atom in circularly polarized microwaves: Chaotic ionization via core scattering. Phys Rev A 1997;55:3730–3745. [4] Celletti, A. The Levi–Civita, KS and Radial–Inversion Regularizing Transformations. LNP 2002;590:25-48. [5] Chenciner, A., Llibre, J. A note on the existence of invariant punctured tori in the planar circular RTBP. Ergod. Th. & Dynam. Sys. 1988;8:63–72. [6] Delgado, J. Transversal ejection-collision orbits in Hill’s problem for C >> 1. Celest Mech Dyn Astron 1988;44(3):299–307. [7] Devaney, R. L. Singularities in Classical Celestial Mechanics. In: Ergodic Theory and Dynamical Systems I, Proceedings Special year, Maryland 1979-80, A. Katok (Ed.), p. 211–333. [8] Dormand, J. R., Prince J. P. A familiy of embedded Runge-Kutta formulae. Journal of Compu- tational and Applied Mathematics, volume 6, n 1, 1980. [9] Hénon, M. Exploration numérique du problème restreint I. Masses égales, Orbites périodiques. Ann Astrophys 1965;28:499–511. [10] Hénon, M. Numerical exploration of the Restricted Problem V. Hill’s case: Periodic orbits and Their Stability. Astron Astrophys 1969;1:223–238. [11] Jorba, A. and Zou, M. A software package for the numerical integration of ODE’s by means of high-order taylor methods. Exp Maths 2005;14:99–117. [12] Lacomba, E. A. and Llibre, J. Transversal Ejection-Collision Orbits for the Restricted Problem and the Hill’s Problem with Applications. J Differ Eq 1988;74:69–85. [13] Llibre, J. On the Restricted Three-Body Problem when the Mass Parameter is Small. Celest Mech Dyn Astron 1982;28:83–105. [14] Llibre, J., Martinez-Alfaro, J. Ejection and collision orbits of the spatial RTBP. Celest Mech Dyn Astron 1985;35:113–128. [15] Llibre, J., Pinyol, C. On the Elliptic Restricted Three-Body Problem. Celest Mech Dyn Astron 1990;48:319–345. [16] McGehee, R. Triple Collision in the Collinear Three-Body Problem. Invent Math 1974;27:191– 227. [17] Ollé, M. To and fro motion for the hydrogen atom in a circularly polarized microwave field. Commun Nonlinear Sci Numer Simulat 2018;54:286-301. 27 M. Ollé, Ó. Rodríguez, J. Soler [18] Ollé, M., Rodríguez, O., Soler, J. Ejection-collision orbits in the RTBP. Commun Nonlinear Sci Numer Simulat 2018;55:298-315. [19] Ollé, M., Rodríguez, O., Soler, J. Regularisation in ejection-collision orbits of the RTBP. In Recent Progress in Pure and Applied Mathematics, RSME Springer Series 2019. [20] Oshima, K. Linking low- to high-energy dynamics of invariant manifolds, transit orbits, and singular collision orbits in the planar circular restricted three-body problem. Cel. Mech. Dynam. Astron. 2019;131:53. [21] Pinyol, C. Ejection-collision orbits with the more massive primary in the planar elliptic restricted three-body problem. Celest Mech Dyn Astron 1995;61:315–331. [22] Stiefel, E. L., Scheifele, G. Linear and regular Celestial Mechanics. Springer-verlag, New York; 1971. [23] Szebehely, V. Theory of orbits. Academy Press, Inc., New York; 1967. 28 Analytical and numerical results on families of n-Ejection-collision orbits in the RTBP Appendix A Computation of V and V −1 This appendix shows the computations needed to obtain V and V −1 . According to (53) we need to ∂ û0 ∂t compute and , which are expressions that depend on s and the initial values u00 . Due to (55) ∂u00 ∂u00 ∂ û0 ∂t and (56) it suffices to compute 0 and . All these derivatives are evaluated on ejection orbits, ∂ û0 ∂ û00√ √ i.e. with initial conditions û00 = u00 = (0, 0, 2 2 cos θ0 , 2 2 sin θ0 ), and they will be denoted with the subscript e after a vertical bar or a square bracket. These expressions are functions of s, θ0 and C. ∂ û0 ∂t Computation of and ∂ û00 e ∂ û00 e " # ∂ û0 4 cos θ0 sin θ0  √  2 cos θ sin θ  √  0 0 • = 1− s cos 2 Cs + sin 2 Cs ∂ û00 C C 3/2 e ∂ û0 4 cos2 θ0  √  2 cos2 θ  √  0 • = s cos 2 Cs − sin 2 Cs ∂v̂00 C C 3/2 e ∂ û0 1  √  • 0 = √ sin 2 Cs ∂ û00 2 C e ∂ û0 • 0 =0 ∂v̂00 e ∂v̂ 0 4 sin2 θ0  √  2 sin2 θ  √  0 • =− s cos 2 Cs + sin 2 Cs ∂ û00 C C 3/2 e " # ∂v̂ 0 4 cos θ0 sin θ0  √  2 cos θ sin θ  √  0 0 • = 1+ s cos 2 Cs − sin 2 Cs ∂v̂00 C C 3/2 e ∂v̂ 0 • 0 =0 ∂ û00 e ∂v̂ 0 1  √  • 0 = √ sin 2 Cs ∂v̂00 2 C e " # ∂ û0 0 √ 8 cos θ0 sin θ0  √  • = −2 C + √ s sin 2 Cs ∂ û00 C e ∂ û0 0 8 cos2 θ0  √  • =− √ s sin 2 Cs ∂v̂00 C e ∂ û0 0  √  • 0 = cos 2 Cs ∂ û00 e 0 ∂ û0 • 0 =0 ∂v̂00 e 29 M. Ollé, Ó. Rodríguez, J. Soler ∂v̂ 0 0 8 sin2 θ0  √  • = √ s sin 2 Cs ∂ û00 C e " # ∂v̂ 0 0 √ 8 cos θ0 sin θ0  √  • =− 2 C+ √ s sin 2 Cs dv̂00 C e 0 ∂v̂ 0 • 0 =0 ∂ û00 e ∂v̂ 0 0  √  • 0 = cos 2 Cs ∂v̂00 e √ √  √  3√2 sin θ  √  2√2 cos θ  √  ∂t 8 2 sin θ0 4 2 sin θ0 0 0 • = s+ s cos 4 Cs − sin 4 Cs + sin2 2 Cs ∂ û00 C 2 C 2 C 5/2 C e √ √  √  3√2 cos θ  √  2√2 sin θ  √  ∂t 8 2 cos θ0 4 2 cos θ0 0 0 2 • =− s− s cos 4 Cs + sin 4 Cs + sin 2 Cs ∂v̂00 C2 C2 C 5/2 C e √ √  √  ∂t 2 2 cos θ0 2 cos θ0 • 0 = s− sin 4 Cs ∂ û00 C 2C 3/2 e √ √  √  ∂t 2 2 sin θ0 2 sin θ0 • 0 = s− sin 4 Cs ∂v̂00 C 2C 3/2 e Computation of V = (Vij ) , 1 ≤ i, j ≤ 4.   ∂ û0 v̂ 0 ∂t • V11 = + ∂u00 2 ∂u00 e  √  4  √  2  √  = cos 2 Cs − sin θ0 cos θ0 cos 2 Cs s + 3/2 sin θ0 cos θ0 sin3 2 Cs C C 4 h  √  i  √  6  √   √  + 5/2 sin2 θ0 2 cos2 2 Cs + 1 sin 2 Cs s − 3 sin2 θ0 sin2 2 Cs cos 2 Cs C C  0  ∂ û v̂ 0 ∂t • V12 = + ∂v00 2 ∂v0 e 4  √  2 h  √ i  √  = cos2 θ0 cos 2 Cs s − 3/2 1 − sin2 θ0 sin2 2 Cs sin 2 Cs C C 4 h  √ i  √  6  √   √  − 5/2 cos θ0 sin θ0 1 + 2 cos2 2 Cs sin 2 Cs s + 3 cos θ0 sin θ0 sin2 2 Cs cos 2 Cs C C  0  ∂ û v̂ 0 ∂t • V13 = 0 + 0 ∂u00 2 ∂u00 e 1  √  2  √  1  √   √  = √ sin 2 Cs + 3/2 cos θ0 sin θ0 sin 2 Cs s − 2 cos θ0 sin θ0 cos 2 Cs sin2 2 Cs 2 C C C  0  ∂ û v̂ 0 ∂t • V14 = 0 + 0 ∂v00 2 ∂v00 e 2  √  1  √   √  = 3/2 sin2 θ0 sin 2 Cs s − 2 sin2 θ0 cos 2 Cs sin2 2 Cs C C 30 Analytical and numerical results on families of n-Ejection-collision orbits in the RTBP   ∂v̂ 0 û0 ∂t • V21 = − ∂u00 2 ∂u00 e 4  √  2 h  √ i  √  = − sin2 θ0 cos 2 Cs s + 3/2 1 − cos2 θ0 sin2 2 Cs sin 2 Cs C C 4 h  √ i  √  6  √   √  − 5/2 cos θ0 sin θ0 1 + 2 cos2 2 Cs sin 2 Cs s + 3 cos θ0 sin θ0 sin2 2 Cs cos 2 Cs C C  0  ∂v̂ û0 ∂t • V22 = 0 − ∂v0 2 ∂v00 e  √  4  √  2  √  = cos 2 Cs + sin θ0 cos θ0 cos 2 Cs s − 3/2 sin θ0 cos θ0 sin3 2 Cs C C 4 h  √  i  √  6  √   √  + 5/2 cos θ0 2 cos 2 Cs + 1 sin 2 Cs s − 3 cos2 θ0 sin2 2 Cs cos 2 Cs 2 2 C C  0  ∂v̂ û0 ∂t • V23 = 0 − 0 ∂u00 2 ∂u00 e 2  √  1  √   √  = − 3/2 cos2 θ0 sin 2 Cs s + 2 cos2 θ0 cos 2 Cs sin2 2 Cs C C  0  ∂v̂ û0 ∂t • V24 = 00 − 0 ∂v0 2 ∂v00 e 1  √  2  √  1  √   √  = √ sin 2 Cs − 3/2 cos θ0 sin θ0 sin 2 Cs s + 2 cos θ0 sin θ0 cos 2 Cs sin2 2 Cs 2 C C C " # ∂ û0 0  2  0 0 0 2 2 v̂ 0 − 2(û0 + v̂ 0 )û0 ∂t 0 02 ∂v̂ 0 0 ∂ û • V31 = + 2 û + 3v̂ + 4û v̂ + ∂u00 ∂u00 ∂u00 2 ∂u00 e √  √  8  √  = −2 C sin 2 Cs + √ cos θ0 sin θ0 sin 2 Cs s C 4 h  √ i  √  + cos θ0 sin θ0 2 − 3 cos2 2 Cs cos 2 Cs C 8 h  √ i  √  − 2 sin2 θ0 5 − 8 cos2 2 Cs cos 2 Cs s C 4 h   √   √ i  √  + 5/2 6 sin2 θ0 − 9 − 11 cos2 θ0 cos2 2 Cs − 2 cos2 θ0 cos4 2 Cs sin 2 Cs C 16 h  √ i  √  − 7/2 cos θ0 sin θ0 1 + 2 cos2 2 Cs sin3 2 Cs s C 24  √   √  + 4 cos θ0 sin θ0 cos 2 Cs sin4 2 Cs C " # ∂ û0 0  2  0 0 0 2 2 v̂ 0 − 2(û0 + v̂ 0 )û0 ∂t 0 02 ∂v̂ 0 0 ∂ û • V32 = + 2 û + 3v̂ + 4û v̂ + ∂v00 ∂v00 ∂v00 2 ∂v00 e 8  √  2 = − √ cos θ0 sin 2 Cs s C 4 h   √ i  √  + 3 − 2 cos2 θ0 − 4 − 3 cos2 θ0 cos2 2 Cs cos 2 Cs C 8 h  √ i  √  + 2 cos θ0 sin θ0 5 − 8 cos2 2 Cs cos 2 Cs s C 4 h  √   √ i  √  − 5/2 cos θ0 sin θ0 6 − 11 cos2 2 Cs + 2 cos4 2 Cs sin 2 Cs C 16 h  √ i  √  + 7/2 cos2 θ0 1 + 2 cos2 2 Cs sin3 2 Cs s C 24  √   √  − 4 cos2 θ0 cos 2 Cs sin4 2 Cs C 31 M. Ollé, Ó. Rodríguez, J. Soler " # ∂ û0 0  2  ∂v̂ 0 0 0 2 2 v̂ 0 − 2(û0 + v̂ 0 )û0 ∂t 0 02 0 0 ∂ û • V33 = 0 + 2 û + 3v̂ 0 + 4û v̂ 0 + 0 ∂u00 ∂u00 ∂u00 2 ∂u00 e  √  4  √  = cos 2 Cs + cos θ0 sin θ0 cos 2 Cs s C 2 h  √ i  √  + 3/2 cos θ0 sin θ0 2 − 3 cos2 2 Cs sin 2 Cs C 8  √  − 5/2 cos2 θ0 sin3 2 Cs s C 4  √   √  + 3 cos2 θ0 cos 2 Cs sin4 2 Cs C " # ∂ û0 0  2  ∂v̂ 0 0 0 2 2 v̂ 0 − 2(û0 + v̂ 0 )û0 ∂t 0 02 0 0 ∂ û • V34 = 0 + 2 û + 3v̂ 0 + 4û v̂ 0 + 0 ∂v00 ∂v00 ∂v00 2 ∂v00 e 4  √  = sin2 θ0 cos 2 Cs s C 2 h   √ i  √  − 3/2 3 − 2 sin2 θ0 − 4 − 3 sin2 θ0 cos2 2 Cs sin 2 Cs C 8  √  − 5/2 cos θ0 sin θ0 sin3 2 Cs s C 4  √   √  + 3 cos θ0 sin θ0 cos 2 Cs sin4 2 Cs C " # ∂v̂ 0 0  2  0 0 0 2 2 û0 + 2(û0 + v̂ 0 )v̂ 0 ∂t 0 02 ∂ û 0 0 ∂v̂ • V41 = − 2 3û + v̂ − 4û v̂ − ∂u00 ∂u00 ∂u00 2 ∂u00 e 8  √  2 = √ sin θ0 sin 2 Cs s C 4 h   √ i  √  − 3 − 2 sin2 θ0 − 4 − 3 sin2 θ0 cos2 2 Cs cos 2 Cs C 8 h  √ i  √  + 2 cos θ0 sin θ0 5 − 8 cos2 2 Cs cos 2 Cs s C 4 h  √   √ i  √  − 5/2 cos θ0 sin θ0 6 − 11 cos2 2 Cs + 2 cos4 2 Cs sin 2 Cs C 16 h  √ i  √  − 7/2 sin2 θ0 1 + 2 cos2 2 Cs sin3 2 Cs s C 24  √   √  + 4 sin2 θ0 cos 2 Cs sin4 2 Cs C " # ∂v̂ 0 0  2  0 0 0 2 2 û0 + 2(û0 + v̂ 0 )v̂ 0 ∂t 0 02 ∂ û 0 0 ∂v̂ • V42 = − 2 3û + v̂ − 4û v̂ − ∂v00 ∂v00 ∂v00 2 ∂v00 e √  √  8  √  = −2 C sin 2 Cs − √ cos θ0 sin θ0 sin 2 Cs s C 4 h  √ i  √  − cos θ0 sin θ0 2 − 3 cos2 2 Cs cos 2 Cs C 8 h  √ i  √  − 2 cos2 θ0 5 − 8 cos2 2 Cs cos 2 Cs s C 4 h   √   √ i  √  + 5/2 6 cos2 θ0 − 9 − 11 sin2 θ0 cos2 2 Cs − 2 sin2 θ0 cos4 2 Cs sin 2 Cs C 16 h  √ i  √  + 7/2 cos θ0 sin θ0 1 + 2 cos2 2 Cs sin3 2 Cs s C 24  √   √  − 4 cos θ0 sin θ0 cos 2 Cs sin4 2 Cs C 32 Analytical and numerical results on families of n-Ejection-collision orbits in the RTBP " # ∂v̂ 0 0  2  ∂ û0 0 0 2 2 û0 + 2(û0 + v̂ 0 )v̂ 0 ∂t 0 02 0 0 ∂v̂ • V43 = 0 − 2 3û + v̂ 0 − 4û v̂ 0 − 0 ∂u00 ∂u00 ∂u00 2 ∂u00 e 4  √  = − cos2 θ0 cos 2 Cs s C 2 h   √ i  √  − 3/2 3 − 2 cos2 θ0 − 4 − 3 cos2 θ0 cos2 2 Cs sin 2 Cs C 8  √  − 5/2 cos θ0 sin θ0 sin3 2 Cs s C 4  √   √  + 3 cos θ0 sin θ0 cos 2 Cs sin4 2 Cs C " # ∂v̂ 0 0  2  ∂ û0 0 0 2 2 û0 + 2(û0 + v̂ 0 )v̂ 0 ∂t 0 02 0 0 ∂v̂ • V44 = 0 − 2 3û + v̂ 0 − 4û v̂ 0 − 0 ∂v00 ∂v00 ∂v00 2 ∂v00 e  √  4  √  = cos 2 Cs − cos θ0 sin θ0 cos 2 Cs s C 2 h  √ i  √  − 3/2 cos θ0 sin θ0 2 − 3 cos2 2 Cs sin 2 Cs C 8  √  − 5/2 sin2 θ0 sin3 2 Cs s C 4  √   √  + 3 sin2 θ0 cos 2 Cs sin4 2 Cs C  Computation of V −1 = Vij−1 , 1 ≤ i, j ≤ 4. According to the theorem of Liouville, the determinant of the matrix V (s) is identically 1 because the trace of DF0 (u0 ) in (16) is zero.  √  4  √  2 h  √ i  √  −1 • V1,1 = cos 2 Cs + cos θ0 sin θ0 cos 2 Cs s + 3/2 2 − 3 cos2 2 Cs cos θ0 sin θ0 sin 2 Cs C C 24  √  12  √   √  + 5/2 cos θ0 sin 2 Cs s − 3 cos θ0 cos 2 Cs sin4 2 Cs 2 3 2 C C 4  √  2 h  √ i  √  −1 • V1,2 =− cos2 θ0 cos 2 Cs s + 3/2 1 + 2 sin2 θ0 − 3 sin2 θ0 cos2 2 Cs sin 2 Cs C C 24  √  12  √   √  + 5/2 cos θ0 sin θ0 sin 2 Cs s − 3 cos θ0 sin θ0 cos 2 Cs sin4 2 Cs 3 C C 1  √  2  √  1  √   √  −1 • V1,3 = − √ sin 2 Cs − 3/2 cos θ0 sin θ0 sin 2 Cs s + 2 cos θ0 sin θ0 cos 2 Cs sin2 2 Cs 2 C C C 2  √  1  √   √  −1 • V1,4 = cos2 θ0 sin 2 Cs s − 2 cos2 θ0 cos 2 Cs sin2 2 Cs C 3/2 C 4  √  2 h  √ i  √  −1 • V2,1 = sin2 θ0 cos 2 Cs s − 3/2 1 + 2 cos2 θ0 − 3 cos2 θ0 cos2 2 Cs sin 2 Cs C C 24  √  3 + 5/2 cos θ0 sin θ0 sin 2 Cs s C 12  √   √  − 3 cos θ0 sin θ0 cos 2 Cs sin4 2 Cs C 33 M. Ollé, Ó. Rodríguez, J. Soler  √  4  √  −1 • V2,2 = cos 2 Cs − cos θ0 sin θ0 cos 2 Cs s C 2 h 2  √ i  √  − 3/2 2 − 3 cos 2 Cs cos θ0 sin θ0 sin 2 Cs C 24  √  + 5/2 sin2 θ0 sin3 2 Cs s C 12  √   √  − 3 sin2 θ0 cos 2 Cs sin4 2 Cs C 2  √  1  √   √  −1 • V2,3 =− sin2 θ0 sin 2 Cs s + 2 sin2 θ0 cos 2 Cs sin2 2 Cs C 3/2 C 1  √  2  √  1  √   √  −1 • V2,4 = − √ sin 2 Cs + 3/2 cos θ0 sin θ0 sin 2 Cs s − 2 cos θ0 sin θ0 cos 2 Cs sin2 2 Cs 2 C C C √  √  8  √  −1 • V3,1 = 2 C sin 2 Cs − √ cos θ0 sin θ0 sin 2 Cs s C 4 h 2  √ i  √  − 2 − 3 cos 2 Cs cos θ0 sin θ0 cos 2 Cs C 8  √  − 2 sin2 θ0 cos 2 Cs s C 4 h   √   √ i  √  − 5/2 6 cos2 θ0 − 1 + 11 cos2 θ0 cos2 2 Cs + 6 cos2 θ0 cos4 2 Cs sin 2 Cs C 48 h  √ i  √  + 7/2 1 − 2 cos2 2 Cs cos θ0 sin θ0 sin3 2 Cs s C 24  √   √  + 4 cos θ0 sin θ0 cos 2 Cs sin4 2 Cs C 8  √  4 h  √ i  √  −1 • V3,2 = − √ sin2 θ0 sin 2 Cs s − 2 − 3 cos2 2 Cs sin2 θ0 cos 2 Cs C C 8  √  + 2 cos θ0 sin θ0 cos 2 Cs s C 4 h  √   √ i  √  − 5/2 6 − 11 cos2 2 Cs + 6 cos4 2 Cs cos θ0 sin θ0 sin 2 Cs C 48 h  √ i  √  + 7/2 1 − 2 cos2 2 Cs sin2 θ0 sin3 2 Cs s C 24  √   √  + 4 sin2 θ0 cos 2 Cs sin4 2 Cs C  √  4  √  −1 • V3,3 = cos 2 Cs − cos θ0 sin θ0 cos 2 Cs s C 2  √  + 3/2 cos θ0 sin θ0 sin3 2 Cs C 4 h  √ i  √  − 5/2 1 − 2 cos2 2 Cs sin2 θ0 sin 2 Cs s C 2  √   √  − 3 sin2 θ0 cos 2 Cs sin2 2 Cs C 4  √  2  √  −1 • V3,4 =− sin2 θ0 cos 2 Cs s − 3/2 cos2 θ0 ∗ sin3 2 Cs C C 4 h 2  √ i  √  + 5/2 1 − 2 cos 2 Cs cos θ0 sin θ0 sin 2 Cs s C 2  √   √  + 3 cos θ0 sin θ0 cos 2 Cs sin2 2 Cs C 34 Analytical and numerical results on families of n-Ejection-collision orbits in the RTBP 8  √  4 h  √ i  √  −1 • V4,1 = √ cos2 θ0 sin 2 Cs s + 2 − 3 cos2 2 Cs cos2 θ0 cos 2 Cs C C 8  √  + 2 cos θ0 sin θ0 cos 2 Cs s C 4 h  √   √ i  √  − 5/2 6 − 11 cos2 2 Cs + 6 cos4 2 Cs cos θ0 sin θ0 sin 2 Cs C 48 h  √ i  √  − 7/2 1 − 2 cos2 2 Cs cos2 θ0 sin3 2 Cs s C 24  √   √  − 4 cos2 θ0 cos 2 Cs sin4 2 Cs C √  √  8  √  −1 • V4,2 = 2 C sin 2 Cs + √ cos θ0 sin θ0 sin 2 Cs s C 4 h  √ i  √  + 2 − 3 cos2 2 Cs cos θ0 sin θ0 cos 2 Cs C 8  √  − 2 cos2 θ0 cos 2 Cs s C 4 h   √   √ i  √  − 5/2 6 sin2 θ0 − 1 + 11 sin2 θ0 cos2 2 Cs + 6 sin2 θ0 cos4 2 Cs sin 2 Cs C 48 h  √ i  √  − 7/2 1 − 2 cos2 2 Cs cos θ0 sin θ0 sin3 2 Cs s C 24  √   √  − 4 cos θ0 sin θ0 cos 2 Cs sin4 2 Cs C 4  √  2  √  −1 • V4,3 = cos2 θ0 cos 2 Cs s + 3/2 sin2 θ0 sin3 2 Cs C C 4 h  √ i  √  + 5/2 1 − 2 cos2 2 Cs cos θ0 sin θ0 sin 2 Cs s C 2  √   √  + 3 cos θ0 sin θ0 cos 2 Cs sin2 2 Cs C  √  4  √  −1 • V4,4 = cos 2 Cs + cos θ0 sin θ0 cos 2 Cs s C 2  √  − 3/2 cos θ0 sin θ0 sin3 2 Cs C 4 h  √ i  √  − 5/2 1 − 2 cos2 2 Cs cos2 θ0 sin 2 Cs s C 2  √   √  − 3 cos2 θ0 cos 2 Cs sin2 2 Cs C Appendix B Computation of u1 √ In this appendix we use the variable τ = 2 Cs introduced in Section 6 in order to simplify the expressions. With this new variable we have from (50) 2τ − sin(2τ ) t(τ ) = (61) C 3/2 The function F1 = (0, 0, F13 , F14 ) in (13) is given by 35 M. Ollé, Ó. Rodríguez, J. Soler √ √ √ 4 2 32 2 8 2 sin τ U1 (θ0 , τ ) • F13 = √ sin τ U1 (θ0 , τ ) + 3/2 sin3 τ U13 (θ0 , τ ) + √ C C Cr2 ! √ 2 16 2 sin3 τ sin2 τ + 1 U1 (θ0 , τ ) C − C 3/2 r23 √ √ √ 4 2 32 2 3 3 8 2 sin τ U2 (θ0 , τ ) • F14 = √ sin τ U2 (θ0 , τ ) − 3/2 sin τ U2 (θ0 , τ ) + √ C C Cr2 ! √ 2 16 2 sin3 τ sin2 τ − 1 U2 (θ0 , τ ) C − C 3/2 r23 where     1 −2τ + sin(2τ ) 1 −2τ + sin(2τ ) • U1 (θ0 , τ ) = cos θ0 cos − sin θ0 sin C 3/2 2 C 3/2 2     1 −2τ + sin(2τ ) 1 −2τ + sin(2τ ) • U2 (θ0 , τ ) = cos θ0 sin + sin θ0 cos C 3/2 2 C 3/2 2 s   2 −2τ + sin(2τ ) 4 • r2 = 1 + cos 2θ0 + 3/2 sin2 τ + 2 sin4 τ C C C The second step is to compute X −1 F1 and its integral over [0, 2k−1 2 π]. The full function F1 gives rise to expressions which are not easy to be integrated. Expanding F1 in powers of C −1/2 up to order C −11/2 we are able to give an analytic result for C large enough. Considering the expansions X (2τ − sin(2τ )) 2i • cos(t(τ )/2) = (−1)i i=0 4i C 3i (2i)! X (2τ − sin(2τ )) 2i+1 • sin(±t(τ )/2) = ∓ (−1)i i=0 22i+1 C 3(2i+1)/2 (2i + 1)! Z 1 we obtain the following expressions for √ X −1 F1 2 C 36 Analytical and numerical results on families of n-Ejection-collision orbits in the RTBP Z (2k−1)π √ √ 2 1 −1  3 2(2k − 1)π cos θ0 15 2(2k − 1)π cos θ0 (2 cos4 θ0 − 1) √ X F1 1 dτ = − − 0 2 C 4C 3/2 4C 7/2 −9/2 + O(C ) Z (2k−1)π √ √ 2 1  3 2(2k − 1)π sin θ0 15 2(2k − 1)π sin θ0 (2 sin4 θ0 − 1) √ X −1 F1 2 dτ = − − 0 2 C 4C 3/2 4C 7/2 + O(C −9/2 ) Z (2k−1)π √ √ √ 2 1 −1  3 2 cos θ0 3 2π sin θ0 (2k − 1) 8 2 cos θ0 (2 cos4 θ0 − 1) √ X F1 3 dτ = + + 0 2 C C 2C 5/2 C3 √ 16 2 cos θ0 (5 cos6 θ0 − 6 cos2 θ0 + 2) − √ C4 2π sin θ0 (2k − 1)(992 sin4 θ0 − 1376 sin2 θ0 + 431) + √ 8C 9/2 32 2 cos θ0 (14 cos θ0 − 30 cos4 θ0 + 20 cos2 θ0 − 3) 8 + + O(C −11/2 ) C5 Z (2k−1)π √ √ √ 2 1 −1  3 2 sin θ0 3 2π cos θ0 (2k − 1) 8 2 sin θ0 (2 sin4 θ0 − 1) √ X F1 4 dτ = − + 0 2 C C 2C 5/2 C3 √ 6 2 16 2 sin θ0 (5 sin θ0 − 6 sin θ0 + 2) + √ C4 2π cos θ0 (2k − 1)(992 cos4 θ0 − 1376 cos2 θ0 + 431) − √ 8C 9/2 32 2 sin θ0 (14 sin θ0 − 30 sin4 θ0 + 20 sin2 θ0 − 3) 8 + + O(C −11/2 ) C5 (62) The final step is to compute u1 (s) from (17). The following formulas are obtained for u1 (s0k ) (actually only the components u1 and v 1 are required) √ √ 1 (−1)k 2 cos θ0 3(−1)k 2 cos θ0 u (s0k ) = √ − 2 C 2C 3/2 k √ √ 3(−1) (2k − 1) 2π sin θ0 3(−1)k (2k − 1) 2π sin θ0 + − √ 4C 2 2C 3 k 2 2 4 (−1) 2 cos θ0 (5π (2k − 1) + 128 cos θ0 − 64) − √ 16C 7/2 (−1) 2 cos θ0 (27(2k − 1)2 π 2 + 1920 cos6 θ0 − 2304 cos2 θ0 + 768) k + + O(C −6 ) 48C 9/2 √ √ (63) 1 0 (−1)k 2 sin θ0 3(−1)k 2 sin θ0 v (sk ) = √ − 2 C 2C 3/2 k √ √ 3(−1) (2k − 1) 2π cos θ0 3(−1)k (2k − 1) 2π cos θ0 − + √ 4C 2 2C 3 k 2 2 4 (−1) 2 sin θ0 (5π (2k − 1) + 128 sin θ0 − 64) − √ 16C 7/2 (−1) 2 sin θ0 (27(2k − 1)2 π 2 + 1920 sin6 θ0 − 2304 sin2 θ0 + 768) k − + O(C −6 ) 48C 9/2 37