Generalized Analytical Results on n-Ejection–Collision Orbits in the RTBP. Analysis of Bifurcations

In the planar circular restricted three-body problem and for any value of the mass parameter μ∈(0,1)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu \in (0,1)$$\end{document} and n≥1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n\ge 1$$\end{document}, we prove the existence of four families of n-ejection–collision (n-EC) orbits, that is, orbits where the particle ejects from a primary, reaches n maxima in the (Euclidean) distance with respect to it and finally collides with the primary. Such EC orbits have a value of the Jacobi constant of the form C=3μ+Ln2/3(1-μ)2/3\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C=3\mu +Ln^{2/3}(1-\mu )^{2/3}$$\end{document}, where L>0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L>0$$\end{document} is big enough but independent of μ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu $$\end{document} and n. In order to prove this optimal result, we consider Levi-Civita’s transformation to regularize the collision with one primary and a perturbative approach using an ad hoc small parameter once a suitable scale in the configuration plane and time has previously been applied. This result improves a previous work where the existence of the n-EC orbits was stated when the mass parameter μ>0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu >0$$\end{document} was small enough. Moreover, for decreasing values of C, there appear some bifurcations which are first numerically investigated and afterward explicit expressions for the approximation of the bifurcation values of C are discussed. Finally, a detailed analysis of the existence of n-EC orbits when μ→1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu \rightarrow 1$$\end{document} is also described. In a natural way, Hill’s problem shows up. For this problem, we prove an analytical result on the existence of four families of n-EC orbits, and numerically, we describe them as well as the appearing bifurcations.


Introduction
This paper studies the existence of ejection-collision orbits in the planar circular Restricted three-body problem (RTBP), which describes the motion of a particle (of neglectible mass) under the attraction of two point massive bodies P 1 and P 2 , called primaries, that describe circular orbits around their common center of mass.Introducing a rotating system of coordinates that rotates with the primaries, and using suitable units of length, time and mass, an autonomous system of four ODE of first order is derived, depending on a unique parameter µ ∈ (0, 1), in such a way that the primaries have masses 1 − µ and µ respectively.Such system of ODE has the well known Jacobi first integral (equal to C along each solution) and is a regular system everywhere except when the particle collides with each of the primaries.
n-ejection-collision orbits (n-EC orbits from now on) are orbits which eject from a primary and its distance to it reaches n relative maxima before colliding with it (see Definition 4.a).Since the n-EC orbits are the main target of this paper and collisions between the particle and one primary lead to singularities in the system of ODE, some kind of regularization, that transforms the original system to a new one which is regular at collisions, is necessary.Among the different possible choices, ranging from local to global regularizations (see [5,7,26,27]) we will use along the paper the (local) Levi-Civita regularization [13], because it is conceptually simple, suitable for our theoretical purposes and efficient for numerical simulations.
The main analytical result of this paper is Theorem 1, where we prove that there exists an L such that for L ≥ L and for any value of µ ∈ (0, 1), n ∈ N and the Jacobi constant C = 3µ + Ln 2/3 (1 − µ) 2/3 , there exist four n-EC orbits, and we characterize them.This improves a recent result (see [20]) where the existence of four n-EC orbits ejecting (and colliding) from the big primary (of mass 1 − µ) is proved but only for small enough µ > 0.
To prove this main result we first consider a weaker version in Theorem 2 where we show that for all n ∈ N, there exists a K(n) such that for K ≥ K(n) and for any value of µ ∈ (0, 1) and C = 3µ+K(1−µ) 2/3 , there exist four n-EC orbits.This weaker version also improves the result of [20] since we cover any value of µ so we can eject from (and collide with) any of the primaries, irrespective of its mass.Another improvement is the proof's approach.In the previous paper a perturbative approach for small enough µ > 0 and big enough C was considered.There, the authors computed the series expansion, with respect to the mass parameter µ, of the ejection (collision) manifold.So the explicit analytical expansion, up to certain order, of this manifold integrated up to a suitable Poincaré section Σ (maximum distance to the ejecting primary) was obtained.For suitable number of crossings with Σ, i for the ejection manifold and j for the collision one (with i + j = n + 1), the resulting two curves C + i and C − j were computed.Achieving such curves required some technicalities, in particular, the computation of terms up to order 9 (at least) in such expansions and the expressions of them in the usual polar coordinates (instead of the initial angle θ 0 ).The application of the Implicit Function Theorem (IFT) to analyze the intersection of both curves gave rise to the existence of four n-EC orbits for any n, C big enough and µ > 0 small enough.
In this paper, the perturbative approach considers a suitable small parameter, related with the inverse of the Jacobi constant, regardless of the value of µ.Moreover, instead of computing the two curves C + i and C − j , we consider the angular momentum at the n-th passage with the minimum distance to the primary (the particle ejected from).We characterize an n-EC orbit by the zero value of that angular momentum.This strategy to use the angular momentum simplifies the computations in three directions: first only expansions up to order 6 are required, second obtaining just one function instead of two different curves, and third the parametrization of the angular momentum directly in terms of θ 0 (thus avoiding the technical issue of the transformation to usual polar coordinates).
The second part of the paper focuses on the bifurcations that may appear when doing the continuation of families of n-EC orbits.It is clear that, given any value of µ > 0, and fixed n, we can continue the four families of n-EC orbits for C big enough, from the IFT.According to previous papers ( [21,20]) we will name such families as α n , β n , δ n and γ n .However, as long as C decreases, the IFT does not apply and bifurcations may appear for suitable values of C. We analyze such bifurcations from the analytical expressions obtained in the series expansions for order higher than 6.A rich variety of bifurcations show up.They are discussed and numerically described.
Precisely the results derived from this numerical exploration provides inspiration to obtain the main result: an explicit expression of the bifurcating value of C as Ĉ = 3µ + Ln 2/3 (1 − µ) 2/3 , i.e. we prove that K(n) = Ln 2/3 .Finally, taking µ → 1 gives rise to the Hill problem.Quite naturally the same kind of proof developed previously applies to the Hill problem.So as a corollary we obtain an analytical result that establishes the existence of four families of n-EC in this problem.Moreover the existence of the successive bifurcations when decreasing C for all n ∈ [1,100] are also numerically discussed.
Concerning previous published results on this subject for the circular planar RTBP, we distinguish between analytical and numerical results.Focusing on the theoretical analysis of n-EC orbits, only the case for n = 1 is considered in Llibre [14], Chenciner and Llibre [4] and Lacomba and Llibre [12].The general case n ≥ 1 is studied in [20], but for small enough values of µ > 0. Regarding a numerical approach, and for n = 1, we mention the papers by Bozis [2], Hénon [8,9], where the authors compute some particular EC orbits that naturally appear when doing the continuation of families of periodic orbits.For the general case n ≥ 1, we mention Ollé, Rodriguez and Soler papers [21,22], where the authors compute and analyze the continuation of families of n-EC orbits for n = 1, ...., 25 and discuss the advantages and disadvantages of Levi-Civita's versus McGehee's [15] regularization.Very recently, Ollé, Rodriguez and Soler [23] analyze the global behavior of the whole set of ejection orbits and the dynamical consequences resulting from the interaction between ejection orbits and the Lyapunov periodic orbit around the collinear equilibrium point L 1 .In particular infinitely many (in a chaotic way) EC orbits show up.
We finally remark that the EC orbits appear quite naturally in astronomical applications.Let us mention that EC orbits allow to explain a mechanism of transfer of mass in binary star systems (see [10,16,25,28]), to describe regions of capture of irregular moons by giant planets ( [1]) or to discuss temporary capture ( [24]).Other applications include the probability of crash motion (see [17,18]) or the role of ejection orbits to explain a mechanism for ionization in atomic problems (see [3,19]).
The paper is organized as follows: In Section 2 we recall some basics of the RTBP, we introduce the Levi-Civita coordinates and the new normalized variables that will become useful to prove the existence of the n-EC orbits for any value of µ > 0. Section 3 recalls the topics described in Section 2 but for the Hill problem.In Section 4 we state the two main theorems, Theorem 1 and Theorem 2, concerning the existence of n-EC orbits in the RTBP.We provide the analytical proof of Theorem 2 in Section 5. Section 6 is devoted to numerically analyse the bifurcations of families of n-EC orbits in the RTBP.In Section 7 we provide the analytical proof of Theorem 1. Section 8 is devoted to the Hill problem.
Finally, we observe that all the numerical computations have been done using double precision and the numerical integration of the systems of ODE rely on an own implemented Runge-Kutta (7)8 integrator with an adaptive step size control described in [6] and a Taylor method implemented on a robust, fast and accurate software package in [11].

The planar RTBP and the Levi Civita regularization
As mentioned in the Introduction, we consider the RTPB.In the rotating (synodical) system, the primaries with mass 1 − µ and µ, µ ∈ (0, 1), have positions P 1 = (µ, 0) and P 2 = (µ − 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 where ˙= d/dt and with r 1 = (x − µ) 2 + y 2 and r 2 = (x − µ + 1) 2 + y 2 .So, the equations become singular when r 1 or r 2 → 0.
The main properties of this system used later on are the following (see [27] for details): 1.There exists a first integral, defined by and known as Jacobi integral.
2. System (1) has the symmetry 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 L i , i = 1, 2, 3, and the triangular ones L i , i = 4, 5. On the plane (x, y), L 1,2,3 are located on the x axis, with x L2 < µ − 1 < x L1 < µ < x L3 and L 4,5 forming an equilateral triangle with the primaries.C Li will stand for the value of C at L i , i = 1, ..., 5.

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 In order to deal with the singularity of the primary P 1 = (µ, 0) (r 1 = 0) we will consider the Levi-Civita regularization (see [27]).The well known transformation of coordinates and time is given by: and we remark that, taking µ ∈ (0, 1) we are regularizing the big primary (if µ ∈ (0, 1/2]) or the small one (if µ ∈ [1/2, 1)).In this new system of coordinates, the solutions of system (1) with Jacobi constant equal to C satisfy: where = d/ds and with The system of ODE is now regular everywhere except at the collision with the primary P 2 (r 2 = 0).
We observe that when studying the system of ODE (6), a value of a Jacobi constant C is fixed.Thus to take an initial condition of this system, we will take (u(0), v(0), u (0), v (0), C(0)).Nevertheless, along the paper, we will actually study system (6) removing the last equation in C, and we will consider the corresponding solution for a fixed C and initial condition simply given by (u(0), v(0), u (0), v (0)).
In this new system of variables, the previous properties of the RTBP are translated as: 1. Jacobi Integral: which is regular at the collision with the primary P 1 .In particular (see [27]), the velocity at the position of the first primary (u = 0, v = 0) satisfies: and therefore the velocities at the collision lie in a circle of radius 8(1 − µ). 2. As the Levi-Civita transformation duplicates the configuration space (see Figure 1) the equations of motion satisfy two symmetries, (9a) as a consequence of the duplication of space and (9b) due to (4): 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.See Figure 1.
4. Similarly, given a value of the Jacobi constant C, the Hill's region in variables (u, v) now becomes In particular we will consider values of the Jacobi constant C ≥ C L1 , the value of the Jacobi constant associated to the equilibrium point L 1 .In this way, it will be enough to regularize only the position of P 1 because the Hill's region associated to these values of C avoids collisions with the second primary (assuming the particle moves in a neighbourhood of P 1 ), see Figure 1).

The Hill problem and the Levi-Civita regularization
The Hill problem is a simplified limiting case of the RTBP that allows to study the vicinity of the small primary when this mass tends to 0 (when mass parameter µ is very small or very close to 1).
We can obtain easily the equation of Hill problem making a translation of the small primary (denoted by P h ) to the origin, and rescaling the coordinates by a factor For our purpose we will consider this second case, so the first step is to introduce new variables (x h , y h ) defined by the relation In this way the expression (2) becomes and taking the limit µ → 1 we obtain the Hill's potential Thus the equations of motion are given by The Hill problem also has some interesting properties for our purposes: 1.The system (13) has a first integral defined by where K is related with the Jacobi integral by: 2. The equations ( 13) not only inherit the symmetry of the problem, that is a symmetry with respect to the x h -axis, but also has an extra one with respect to the y h -axis.In this way the system (13) has the symmetries: 3. The Hill problem only preserves two equilibrium points, which are those that are in the vicinity of the small primary P h .That is L 1 and L 2 if we consider µ → 0 or L 1 and L 3 if µ → 1.
For historical consistency, we will call these equilibrium points L 1 and L 2 , which have positions (±1/3 1/3 , 0) (see Figure 2) and we will denote by K L = 3 4/3 the value of K at L 1 and L 2 .
4. In a similar way, from the first integral and taking into account that 2Ψ(x h , y h ) − K ≥ 0, given a value of K, the motion can only take place in the Hill's region defined by We notice that, similarly as we do in the RTBP, we will consider values of K ≥ K L to guarantee that if the particle starts in a region around P h , it will always remain there.
In order to regularize the Hill problem we only have to consider the Levi-Civita regularization and the system (13) becomes: with Under this transformation the previous properties of the Hill problem are translated as: 1.The first integral of ( 18) is given by which is regular at the collision with P h .In particular the velocity at the position of the primary 2. As the transformation duplicates the configuration space (see Figure 2), the equation ( 18) has an extra symmetry: 3. For the same reason, the equilibrium points are duplicated, and we have L 1 = (±3 −1/6 , 0) and L 2 = (0, ±3 −1/6 ).
4. In a similar way, depending on the value of K we can define the valid region of motion (see Figure 2) in the plane (u h , v h ) as: 4 n-EC orbits in the RTBP and main theorems In this paper we will focus on a specific type of EC orbits, the n-EC orbits, formally defined as Definition 4.a.We call n-ejection-collision orbit of a primary, simply noted by n-EC orbit, to the orbit that the particle describes when ejects from a primary and reaches n times a relative maximum in the distance with respect to this primary before colliding with it.
As we will consider any value of µ ∈ (0, 1) we will study only the n-EC orbits associated to the first primary P 1 .Notice that from relation (8) it is easy to compute the initial conditions of the ejection orbits (and the collision orbits): and we can compute the manifold of the ejection (collision) orbits integrating forward (backward) in time.Observe that in this case it is enough to consider a value of θ 0 ∈ [0, π) due to the duplication of the configuration plane.
Concerning the existence of n-EC orbits, we mentioned above that in [20], the existence of four n-EC orbits ejecting from (and colliding with) the big primary for any n ≥ 1, given C big enough and µ > 0 small enough, was proved.The proof was based on a perturbative approach in µ and assuming that the orbits ejected from the big primary of mass 1 − µ.
The first goal of this paper is to improve this previous result and prove the existence of four n-EC orbits ejecting from (and colliding with) the big or small primary, for any n ≥ 1 given and C big enough.So any value of the mass parameter µ ∈ (0, 1) is possible in this context.
For analytical and numerical purposes, though, we will use a characterization for an EC orbit, based upon the zero value of its angular momentum, defined from now on as M := U V − V U (for some suitable variables (U, V ) to be defined later), at a minimum distance with the primary the particle ejected from (see Lemma 1 below).So in order to obtain an n-EC orbit, for n ≥ 1, µ ∈ (0, 1) and C given, first we will compute the corresponding ejection solution for each initial condition (that is, for each value θ 0 ).Second we will determine the precise time τ * = τ * (θ 0 ) when the particle reaches the n-th minimum in the distance to P 1 .At time τ * we will compute the value of the angular momentum that is, (U V − V U )(τ * ).Varying θ 0 ∈ [0, π) we will obtain the corresponding angular momentum, that will denote by M (n, θ 0 ) = (U V − V U )(τ * ) (overlooking the additional dependence on µ).The zeros of M (n, θ 0 ) = 0 will provide us with the precise values of θ 0 such that the corresponding ejection orbit is precisely an n-EC orbit.Just to show this idea, we plot in Figure 3 left the behaviour of the angular momentum M (1, θ 0 ) for µ = 0.1 and C = 5.In the right figure we plot the corresponding ejection orbits for three chosen values of θ 0 : the red one and green one plotted for a range of time [0, τ * + δ] (a small suitable δ > 0) to see the change of sign in the angular momentum (shown in the zoom area) and the blue one which is a 1-EC orbit.Now, we proceed to state the main result of this paper about the existence, the number and the characteristics of the n-ejection-collision orbits for any value of the mass parameter and n ∈ N, for sufficiently restricted Hill regions (i.e.C big enough).
Theorem 1.There exists an L such that for L ≥ L and for any value of µ ∈ (0, 1), n ∈ N and , 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, one symmetric of the other with respect to the x axis.
The respective families are labelled by α n , γ n , β n and δ n (corresponding to increasing values of θ 0 respectively).In order to prove Theorem 1 we will first state a weaker version of this theorem, Theorem 2, but the proof of this second version will provide light on the approach, mainly a suitable scaling in the configuration variables, time and the Jacobi constant C, used to prove the more optimal result in Theorem 1.
Theorem 2. For all n ∈ N, there exists a K(n) such that for K ≥ K(n) and for any value of µ ∈ (0, 1) and , there exist four n-EC orbits, which can be characterized in the same way as in Theorem 1.
We remark that, in Theorem 2, we have a uniform constant K = K(n) for any value of µ ∈ (0, 1).This implies that when µ → 1 the value of the Jacobi constant (for which Theorem 2 holds) tends to 3, as C L1 does as well.Precisely, and as shown in the proof of Theorem 2, the expansion of C L1 was the inspiration to choose a suitable scaling in the variables, time and C. Finally in Theorem 1 an expression for K(n) as Ln 2/3 is provided.

Proof of Theorem 2
In order to prove Theorem 2, let us fix C ≥ C L1 and consider the following change of variables and time: that corresponds to the change that normalizes the linear term of ( 6) and the initial condition of the ejection orbits.Denoting by ˙= d dτ the new time derivative the system (6) transforms to the following: where It is important to remark that the properties ( 7), ( 8), (10) are preserved (translated to the new variables), and so are the symmetries obtained in the Levi-Civita regularization, i. e.: At this point, the two main ideas to prove the theorem are: (i) a perturvative approach taking δ = 1/ √ C − 3µ as a small parameter, and (ii) the requirement of the angular momentum to be zero at a minimum distance with the primary the particle ejected from.
First of all we observe that the functions 1/R 2 and 1/R 3 2 are are analytic for U , V bounded, 0 ≤ µ ≤ 1, and δ small enough.In fact, the expansions of 1/R 2 and 1/R 3  2 are of the form: where P 2k (U, V ) and Q 2k (U, V ) are polynomials sum of monomials of degree 2k.
So if we expand the system (26) with respect to δ we obtain: which is an analytical system of ODE in δ and P2k−1 (U, V ) and Q2k−1 (U, V ) are polynomials sum of monomials of degree 2k − 1.
Before proceeding it is important to make two observations: 1. We can introduce the parameter ε = (1 − µ) 1/3 δ.So we have, using that δ = 1 C−3µ : 2. We also know that C ≥ C L1 (µ) since otherwise the Hill region of motion allows transits between both primaries and, in this sense, Hill's region is not regular anymore.As it is well known the expansion of C L1 (µ) is (see [27]) therefore, if it is possible we will like to have a uniform parameter K with the same order in (1 − µ).So, introducing the variable K as we have that the previous expression (30) becomes: The change ( 25) using (32) becomes: Note that the value K is related with the Hill constant by (15) when µ tends to 1.
So, in terms of ε = 1/ √ K the system (29) has the following expression: Second let us prove the following characterization for an EC orbit, based upon the zero value of the angular momentum at a minimum distance with the primary.
Lemma 1. Assume C large enough.An ejection orbit is an EC orbit if and only if it satisfies that at a minimum in the distance (with the primary) the angular momentum Proof.The minimum distance condition is given by: and the angular momentum condition M = U V − V U = 0: We will distinguish between two cases: 1. V = 0.Then, from (36): and, by (36) also U = 0.
On the other hand, it is clear that if a collision takes place, i. e. U = V = 0 and U 2 + V 2 = 8(1 − µ), then conditions (35) and (36) are trivially satisfied.Now let us proceed.Using the vectorial notation U = (U, V, U , V ) T , the second order system of ODE (34) can be written as where We remark that G 0 and G 3 are the only functions that depend on U and V , the remaining ones depending only on U and V .Moreover we observe that the expressions appearing in the expansions are polynomials.Both properties allow to significantly simplify the computations.
The next natural step consists on obtaining a solution U = U (τ ) as a series expansion in ε: As a usual procedure to obtain the functions U j , we plug U in system (37), and comparing the powers in ε, we obtain a system of ODE for U j .

Computation of the functions U j
Now we proceed to compute the explicit expressions for U j (τ ) = (U j (τ ), V j (τ ), Uj (τ ), Vj (τ )), for any j.Actually we will show that, in order to prove Theorem 2, we only need to find explicitly the functions U j up to order j = 6.
From Definition 4.a and the scaling (33), any ejection orbit U (τ ) = U (τ, θ 0 ), has the initial condition so we have U 0 (0) = (0, 0, cos θ 0 , sin θ 0 ), Solution for ε = 0: We must solve the linear system: which is a harmonic oscillator, with initial condition (40).Then the ejection orbit U 0 is given by Solution for ε = 0: In order to find the functions U j , we must solve the successive resulting ODE when substituting U by the series expansion in (34) up to the desired order.
We observe that, for j ≥ 1, the linear non homogeneous system of ODE to be solved is where the homogeneous system is always the same but the independent term changes and increases in complexity with j.
Since a fundamental matrix for the homogeneous system (the first order variational equations) is given by and the initial conditions are U j (0) = 0 for j ≥ 1, we obtain the following well known formula The corresponding explicit expressions are the following: ).Once we have the ejection solution up to order j = 6, the next step consists of computing the n-th minimum in the distance to the primary (located at the origin) the particle ejected from as a function of the initial θ 0 .Equivalently we want to compute the n-th minimum of the function U 2 + V 2 (τ ).This requires to compute the precise time denoted by τ * , needed to reach the n-th minimum in distance.We apply the Implicit Function Theorem to the function (U U + V V )(τ * ) = 0 in order to obtain an expansion series in ε, i.e.: We can easily compute τ * 0 , since we have a harmonic oscillator: Writing the function (U U + V V )(τ ) as an expansion series in ε and collecting terms of the same order, we can successively find the terms τ * i (up to order 6, higher order terms in Appendix A): with τ i (n, θ 0 ) = 0 for i = 1, 2, 3, 4, 5(, 7).

Now we are ready to compute the angular momentum
or in short, since we look for the zeros of M (n, θ 0 ) = 0, we write, dividing the previous equation by µε 6 , Now we apply the Implicit Function Theorem and for ε > 0 small enough we obtain that the equation (49) has four and only four roots in [0, π) given by regardless of the value of the parameter µ.It is clear from (49) that the roots θ 0 are simple.
So we have proved that there exist four n-EC orbits.Moreover, applying the symmetries of the system we can conclude that those EC orbits with an intersection angle with m = 0, 2 correspond to symmetric n-EC orbits (in the sense that the (x, y) projection is symmetric with respect to the x axis).Such EC orbits belong to families γ n and α n respectively when varying ε (or, equivalently, K and therefore C in (32)).Those EC orbits with an intersection angle with m = 1, 3 correspond to symmetric n-EC orbits (in the sense that the (x, y) projection is symmetric one with respect to the other one).Such EC orbits belong to families δ n and β n respectively.
This finishes the proof of Theorem 2.
In order to illustrate the results of Theorem 2, in Figure 4

Analysis of Bifurcations
So far we have applied the Implicit Function Theorem to infer the existence of four and only four n-EC orbits, for any value of µ and C = 3µ + K(1 − µ) 1/3 (see (32)) with K big enough, that is ε = 1/ √ K small enough.In this procedure the minimum order required in the ε expansions for both the functions U j and τ * j was order 6.Of course, when ε becomes bigger, the Implicit Function Theorem may not be applied anymore and bifurcations can appear.This section is focused on such bifurcations.
We will focus on two purposes: on the one hand, the illustration of the appearance and collapsing of bifurcating families of n-EC orbits when doing the continuation of families varying C as parameter; and on the other hand, the behavior of K(n) and its associated value Ĉ(µ, n) = 3µ + K(n)(1 − µ) 1/3 provided by Theorem 2, for any value of µ ∈ (0, 1) and varying n.

Bifurcating families
The first task is to compute the angular momentum M (n, θ 0 ) to higher order.To do so we need higher order terms for both the functions U j and τ * j .We have proceeded as in the previous Section; however, there, expressions up to order 6 were enough.To analyze the bifurcations, we provide their expressions for j up to order 10 in the Appendix.Now we are ready to compute the explicit expression for the angular momentum M (n, θ 0 ) = (U V − V U )(τ * ) up to order 10 which is the following: It is clear that if ε is small enough, the dominant term is ε 6 , and the zeros of M (n, θ 0 ) are related to the term sin(4θ 0 ).Therefore we obtain four n-EC orbits.We will illustrate two different kind of bifurcations that take place when doing the continuation of families of n-EC orbits and that can be explained precisely from the analytical expression of M (n, θ 0 ) to higher order just obtained.
The first kind of bifurcation can be inferred just taking into account the terms of M (n, θ 0 ) up to order 8 in (51).The bifurcation is associated with the term sin(6θ 0 ).See Figure 5 top for µ = 0.1 and n = 2.We can clearly see how increasing ε (decreasing C), the bifurcation takes place.Let us describe the bifurcation close to the 2-EC orbit belonging to family α 2 .See the zoom area in Figure 5 top.Locally, at a neighbourhood of the value of θ 0 of such EC orbit, for some value of C the angular momentum has a unique transversal intersection with the x-axis (that is M (2, θ 0 ) = 0, M (2, θ 0 ) = 0).For C = 3.76 this intersection corresponds to the 2-EC orbit belonging to the family α 2 (see the blue curve).For the bifurcating value C bif = 3.72442505, M (2, θ 0 ) crosses tangently the x axis (see the red curve).For smaller values of C, M (2, θ 0 ) crosses the x axis three times, giving rise to two new bifurcating families of 2-EC orbits (see the green curve) besides family α 2 which persists.The new 2-EC orbits are (obviously due to symmetry (4)) one symmetric with respect to the other.From a global point of view, for a range C < C bif , varying θ 0 ∈ [0, π), M (2, θ 0 ) crosses six times, that is, we obtain six 2-EC orbits, and this is related to the term sin(6θ 0 ), which becomes the dominant term in M (2, θ 0 ).We show these 2-EC orbits in Figure 5 bottom.More precisely, on the three plots, the four 2-EC orbits are shown (in the plane (x, y)) and those 2-EC orbits of family α 2 are plotted in a darker colour.Since the family α 2 persists after the bifurcation, the 2-EC orbits are plotted in the left, middle and right plots.The two new bifurcating 2-EC orbits after the bifurcation are also shown on the right plot.
The second kind of bifurcation can be inferred from the expression of M (n, θ 0 ) up to order 10 given in (51).The bifurcation is associated with the term sin(8θ 0 ).See Figure 6 top for µ = 0.1 and n = 3.
We can clearly see how increasing ε (decreasing C), the angular momentum M (3, θ 0 ) typically crosses four times the x-axis (for θ 0 ∈ [0, π)), as expected (see the blue curve in the top figure).However at some bifurcating value C bif there appear two tangencies (say from nowhere, see the red curve in the zoom area in Figure 6 top); each tangency gives rise to two families when doing the continuation of families decreasing C. See the green curve in the zoom area in Figure 6 top.So from a global point of view, for a range of C < C bif and θ 0 ∈ [0, π), the angular momentum M (3, θ 0 ) = 0 crosses eight times the x-axis, giving rise to eight 3-EC orbits related to the term sin(8θ 0 ).We show these 3-EC orbits in Figure 6 bottom.More specifically, on the three plots, the four 3-EC orbits are shown (in the plane (x, y)) and those 3-EC orbits of family α 3 are plotted in a darker colour.The two 3-EC orbits that appear due to the tangency of M (3, θ 0 ) with the x-axis are also plotted on the middle plot.Moreover, the four new bifurcating 3-EC orbits after the bifurcation are also shown on the right plot.A continuous and discontinuous line with the same colour correspond to EC orbits that are symmetric one with respect to the other one.In fact, due to the symmetry of the problem, we might only consider the two intersection points (those on the left hand side or on the right one of the value of θ 0 in α 3 ), and the other two intersection points would be obtained by symmetry.
So far we have described two specific kinds of bifurcations that take place for n = 2 and n = 3, for µ = 0.1.But from the expression of the angular momentum (51) and the previous discussion, we can foresee a great and rich variety of bifurcations.To have a global and exhaustive insight, we have done massive numerical simulations in the following sense: we have fixed a value of µ, and, for a range of values of (ii) In the first row, right plot, and C close to 3.7, we see the first kind of bifurcation described above for n = 2.In the second row, left plot, and C close to 3.9, we recognise the second kind of bifurcation described above for n = 3. (iii) It is clear from such diagrams that, when we decrease C and increase the value of n, several phenomena of collapse of families and bifurcation of new families are more visible.See for example the third row plots, when decreasing C, for θ 0 < π/2, the collapse of two families, and the appearance of two new ones for θ 0 > π/2.Even richer are the diagrams on the last row of the figure.(iv) We have also plotted on each bifurcation diagram the value of the first bifurcation value of C (decreasing C), which is precisely the value Ĉ(µ, n), for µ = 0.1 mentioned in Theorem 2. We notice in the plot how this value Ĉ(µ, n) increases when n increases.
When we take a bigger value of µ, for example, µ = 0.8, we obtain Figure 8. Comparing the plots obtained with those of Figure 7, we observe two effects: the value of Ĉ(µ, n) is smaller, for the same value of n, and moreover, for n = 2, 3, 4, a value of Ĉ(µ, n) really smaller than C L1 is required (compare the four first plots in Figures 7 and 8).For bigger values of µ and for the same value of C ≥ C L1 , the the Hill region gets really smaller, when increasing µ, so quite naturally, the probability of bifurcations decreases.On the other hand, taking C < C L1 represents an enlarging of the Hill's region and therefore a more powerful influence of the big primary, so an easier scenario to have bifurcations.
Remark.We notice that Lemma 1 provides a characterization of an EC orbit if C is large enough.
Along the numerical simulations done, where the values of C are not so large, we have also used the same characterization, but additionally checking that when M = 0 at a minimum distance, then U = V = 0.

Behaviour of Ĉ(µ, n)
As a final goal, we want to describe (numerically) the behaviour of Ĉ(µ, n) = 3µ + K(1 − µ) 2/3 for any value of µ ∈ (0, 1) and n.More precisely, for each value of µ and n, and C big enough, Theorem 1 claims that there exist exactly four families of n-EC orbits.As discussed in the previous subsection, when decreasing C bifurcations appear in a natural way.So for fixed µ and n, the first value of C (decreasing C) such that there appear more than four n-EC orbits is precisely the value Ĉ(µ, n) formulated in Theorem 1.In the previous subsection, we have computed the value Ĉ(µ, n), just for µ = 0.1 and n = 1, ..., 8. Our purpose now is to compute Ĉ(µ, n) for any µ ∈ (0, 1) and n.We will always assume that any value of C considered satisfies C ≥ C L1 (recall the Hill regions in Figure 1, there is no possible connection between P 1 and P 2 , and therefore the dynamics around each primary is the simplest possible).
The strategy to compute numerically Ĉ, for a fixed µ ∈ (0, 1) and given n, is the following: we take the interval I = [C L1 , C b ] of values of C, and for each C ∈ I, (starting at C b ) we vary θ 0 ∈ [0, π) (that defines the initial conditions of an ejection orbit in synodical Levi-Civita variables) and find the four specific values of θ 0 (such that M (n, θ 0 ) = 0) corresponding to the expected four n-EC orbits.So we have four n-EC orbits for that value of C and decreasing C we obtain four families of n-EC orbits.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 new families of n-EC orbits.That is precisely the specific value of Ĉ.
In Figure 9 we show the results obtained for µ ∈ (0, 1) and n = 2, ..., 10.Also the curve (µ, C L1 ) has been plotted (in black).Recall that, as mentioned above, we are focussed on values of C ≥ C L1 .We remark that for n = 1, the value Ĉ(µ, 1) is less than C L1 and therefore is not considered.Moreover for the specific values of µ = 0.1 and µ = 0.8 we recover the indicated values in Figures 7 and  8 respectively.From Figure 9 it is clear that the value of Ĉ(µ, n) increases when n increases.This means that for higher values of n, that is longer time spans integrations, the effect of the other primary is more visible.
We remark that the shape of the curves in Figure 9 provide a hint about the dependence of K(n) in the expression Ĉ(µ, n) = 3µ + K(n)(1 − µ) 2/3 obtained in Theorem 2, more specifically K(n) = Ln 2/3 .We will prove rigorously this dependence in Theorem 1 in the next Section.

Proof of Theorem 1
The proof of Theorem 1 is also based on a perturvative approach.let us introduce a new parameter L defined as K = Ln 2/3 in (32).In this way, we perform the change (33) and a new time T = τ /n: where we introduce the functions in the new time: The system (6) becomes, denoting . Let us introduce the parameter ξ = 1/ √ L, in this way the system (53) becomes with Let us introduce the vectorial notation U = (U, V, U, V).The system (54) can be written as where Remark.Note that F 1 (U ) only depends on the position variables, At this point, our next goal is to find the solution as U = U 0 + U 1 where Notice that U 0 is the solution of the 2-body problem (µ = 0) in synodical (rotating) Levi-Civita coordinates.That is, we consider system (56) as a perturbation of the 2-body problem (58a) where the perturbation parameter is ξ which will be small enough and for any value of µ ∈ (0, 1).Roughly speaking, for big values of the Jacobi constant the problem is close the two body problem of the mass-less body and the collision primary, regardless the value of the mass parameter µ.
Note that we are interested only in the ejection orbits U e = U e 0 + U e 1 and the initial conditions of these orbits are given by U e 0 (0) = (0, 0, n cos θ 0 , n sin θ 0 ) and U e 1 (0 To prove the theorem we will use the same strategy of computing the angular momentum M(n, θ 0 ) at the n-th minimum of the distance to the origin and find the values of θ 0 such that M(n, θ 0 ) = 0.
Thus, we will compute U e and the time needed to reach n-th minimum solving U e Ue + V e Ve (θ 0 , T * ) = 0.The last step will be to calculate M(n, θ 0 ).
In order to obtain the solution of this system, we first consider the 2-body problem in sidereal coordinates: and the change of time being ( Ū0 ( T ), V0 ( T ) the associated solutions.
We recall that for the two body problem the angular momentum is constant, consequently and therefore the solution of ( 61) is given by where ω = n 2 − 4( Ū0 V0 − V0 U0 )(0)ξ 3 .Moreover the value t( T ) is simply obtained from ( 62) and (64): Now we apply the rotation transformation to (64) to obtain the solution of system in synodical coordinates (60), Notice that the relation between the sidereal initial conditions and the synodical ones (U 0 V 0 , U0 V0 )(0), are obtained simply from (66) putting T = 0 Since we are interested in the particular case of ejection orbits, which have as their initial condition the corresponding ejection solution is given by U e 0 = (U e 0 , V e 0 , Ue 0 , Ve 0 ), where: with If we denote by T * 0 the time needed by U e 0 ( T ) to reach the n-th minimum distance to the origin, it is clear from (69) that The perturbed system In order to solve the perturbed problem (i.e.µ = 0) we rewrite system (58b) as where U 0 = U e 0 is the ejection solution (69) of the two body problem and Note that the ejection solution U e 1 has zero initial condition and therefore is the solution of the implicit equation where we define and X( T ) is the fundamental matrix of the linear system: We will apply a Fixed Point Theorem to prove the existence of the solution U e 1 .Thus we consider the space χ = {U : [0, T ] −→ R 4 , U continuous}, for a given T , for example T = 2π.
For a given function U = (U, V, U, V) ∈ χ we consider the norm: With this norm χ is a Banach space.
As usual, given an R > 0, we define the ball B R (0) ⊂ χ as the functions U ∈ χ such that U ≤ R.
Next lemmas show that the required hypotheses for the Fixed Point Theorem to be applied are satisfied.

Proof. See Appendix B.1.
Lemma 3.There exist 0 < ξ 1 ≤ ξ 0 and a constant At this point we select ξ 1 s.t.M 2 µξ 6 1 < 1/2, so we have the following result Lemma 4.Under the same hypotheses of Lemma 3 if we reduce ξ 1 such that M 2 µξ 6 1 < 1/2, one has that the operator H : B R (0) → B R (0) and it is a contraction and therefore there exists a unique U e 1 ∈ B R (0) which is solution of the equation (74) in χ.
Proof.If U ∈ B R (0), then: and we already know by Lemma 3 that H is Lipschitz with Lipschitz constant M 2 µξ 6 < 1/2.
By the Fixed Point Theorem there exists a unique U e 1 ∈ B R (0) which is solution of the equation (74).
Observe that once we know the existence and bounds of the function U e 1 , its smoothness is a consequence of being solution of a smooth differential equation.
The results of the previous lemmas give us the following properties: Writing these inequalities in components, and using the definition of the norm (77), we have • Ue where Lemma 5.With the same hypotheses of Lemma 4, the value of H{0}( T * 0 ) = H{0}(π) is given by where U e 6 ( T * 0 ) = (U e 6 , V e 6 , Ue 6 , Ve 6 )( T * 0 ) are the coefficients of U e 1 of order 6 in ξ evaluated at T * 0 .They are given by: Proof.See Appendix B.3.
With this notation, we have The time needed to reach the n minimum in the distance with the first primary can be obtained from the following Lemma: Lemma 6.With the same hypotheses of Lemma 4, the time T * needed for the ejection solution U e to reach the n minimum in the distance with the first primary is given by T * = T * 0 + T * 1 , where T * 0 = π and: Proof.In order to compute the n minimum in the distance with the first primary we have to solve and therefore we have Finally the angular momentum at T * is given by Lemma 7.With the same hypotheses of Lemma 4, the angular momentum of the ejection solution U e at time T * is given by: Proof.In this way, applying the Implicit Function Theorem, we have that for ξ ≥ 0 small enough we obtain that M(n, θ 0 ) has four and only four roots in [0, π) given by regardless of the values of the mass parameter µ and n.We can characterize this n-EC orbits in the same way as in Theorem 2.
This concludes the proof of Theorem 1.

Results for the Hill problem
As we have seen in Section 3, Hill problem is a limit case of RTBP.In this way, the results obtained in the previous sections can easily be extrapolated to the case of Hill problem.In particular, if in (18)  Furthermore, thanks to the fact that the polynomials P2k and Q2k disappear, it is not necessary to consider an expansion in terms of ε = 1/ √ K and it can be considered directly an expansion on = 1/K 3/2 .
Similarly, if we introduce K = Ln 2/3 , that is, we consider the change we obtain the same system of equations as (54) putting µ = 1 and considering ξ = 1/ √ L. In this way we can obtain the following corollary of Theorem 1: Corollary 8.2.There exists an L such that for L ≥ L and for any value of n ∈ N and K = Ln 2/3 , there exist exactly four n-EC orbits, which can be characterized in the same way as the previous corollary.
In this way, if we do the numerical exploration to compute the n-EC orbits that exist for values of K ≥ K L (see Figure 11) we see that, as expected by the Corollary 8.2, the value of K grows with n.
Before going into more detail on the value of K let us make a few comments about Figure 11.It is important to note that thanks to the extra symmetry we could only study the ejection orbits with θ 0 ∈ [0, π/2), but in order to visualize the evolution of the n-EC orbits we will consider the interval θ 0 ∈ [0, π) in Figure 11.In this figure we observe how at least the first new families of n-EC orbits that appear are born from two of the original families (α n and γ n , or β n and δ n ) when the angle of ejection θ 0 is 0 and π/2 respectively (i.e.ϑ 0 = 0, π) and collapse into the two other original families when the value of θ 0 is π/4 and 3π/4 (i.e.ϑ 0 = π/2, 3π/2) (see for example Figure 12).These respective values are very particular, since when these bifurcations take place we have that the n-EC orbits are periodic or are part of a periodic EC orbit.In particular we have: • If the θ 0 of β n is 0 or π/2 (therefore θ 0 of δ n is π/2 or 0) then we have periodic EC orbit formed by β n and δ n (see Figure 12 left).Analogously, if the θ 0 of α n is π/4 or 3π/4 (therefore θ 0 of γ n is 3π/4 or π/4) then we have periodic EC orbit formed by α n and γ n (see Figure 12 right).
We have computed the value K(n) for n = 1, . . ., 100.See Figure 14).It is important to remark that the numerical value of K(n) obtained has the same shape as the analytical bound of the Corollary 8.2.In particular, if we draw the curve Ln 2/3 with L = 2 2/3 we can see how it practically matches the value of the numerical bound obtained for K (see Figure 14).
To conclude, we have seen how not only does the value of K(n) follow the curve Ln 2/3 with L = 2 2/3 , but also the successive bifurcations (the values of K where appear new EC orbits) are closely related to the curves Ln 2/3 with L = (2/p) 2/3 being p a natural number.In particular, in Figure 15 we can see how the value of the successive bifurcations coincides with the curves Ln 2/3 with L = (2/p) 2/3 and p = 1, ..., 10.Then, writing the function U U + V V as an expansion series in ε and collecting terms of the same order, we can successively find the terms τ * i of order i = 7, ..., 10 from (U U + V V )(τ * ) = 0: Now we are ready to compute the explicit expression for the angular momentum M (n, θ 0 ) = (U V − V U )(τ * ) up to order 10 which is the following: which is precisely (51).
During this section we will use M to denote any constant which appears in the bounds and is independent of ξ, µ and n ∈ N.
Lemma 8.The fundamental matrix X for system (76) can be expressed as , and its inverse matrix as Proof.Consider the general solution U 0 of system (58a) given by ( 64), ( 65) and (66), we can express the fundamental matrix of the system and A is the matrix with rows where Ū 0 is given by (64).
Note that we are interested in solving equations (76), which correspond to the ejection orbits U e 0 , that have initial conditions (0, 0, n cos θ 0 , n sin θ 0 ).So we must compute the fundamental matrix X with U 0 ( T ) = U e 0 ( T ).We denote by A e and R e the corresponding matrices.Recall also that the expression of t is given by (70) and the explicit elements of A e are provided in Appendix C.
In this way, we can express A e and R e as and therefore, The expression for X −1 ( T ) can be found in a similar way.

B.1 Proof of Lemma 2
From (73) and (75) we have so, the first step is to bound the components of F 1 (U e 0 ) (see ( 57)).
Concerning the expansions involving r 2 in (55), we have where the symbol O refers to terms bounded for bounded U and any µ ∈ (0, 1) and n ∈ N. Thus, we obtain and the constant M is independent of µ and n.

By Lemma 8 we can bound |X| ≤ M and |X
In this way, we have And, therefore, as we have taken T = 2π: Finally, multiplying by µX we have and using the definition of the norm in (77) and renaming M 1 = M we obtain the desired result:
Lemma 9. Take U ⊕ , U ∈ B R (0).Then for 0 < ξ small enough we have that: Proof.First, we observe that G 0 (U ) = (G 1 0 , G 2 0 , G 3 0 , G 4 0 )(U ) = (0, 0, G 3 0 , G 4 0 )(U ).Therefore we will consider the last two components.We will do the computations for G 3 0 , the ones for G 4 0 are analogous.Using the Mean Value Theorem we have: (100) Taking into account the integral expression in (98) we obtain We get a similar bound for the fourth components and using that the first and the second are identically zero and the definition of the norm we get the result of the lemma.
The next goal is to bound G 1 (U ⊕ ) − G 1 (U ).To do so, we apply the same trick: Lemma 10.Given U ⊕ , U ∈ B R (0).Then for ξ > 0 small enough we have that Proof.Using again the Main Value Theorem we obtain: Let us recall that DG 1 (U ) = µDF 1 (U e 0 + U , V e 0 + V ), and F 1 is given in (57).Proceeding similarly as to bound F 1 we can differentiate (88) to obtain: And using the integral equation ( 101) and the fact that the first two rows of the previous matrix are zero we get: Now, using the definition of the norm we get the result.
From the results of lemmas 9 and 10 we have: Now, proceeding as we did in the proof of Lemma 2, we multiply G(U ⊕ ) − G(U ) by X −1 , integrate for a finite time and multiply the resulting expression by X.We use that |X| ≤ M and |X −1 | ≤ M where M is given in (91) and proceed to bound the expression which gives H{U ⊕ } − H{U } as we did for H{0} in (92), (93), (94), (95), to obtain This finishes the proof of Lemma 3.

B.3 Proof of Lemma 5
In order to compute H{0}( T * 0 ) let us recall its expression:

C Value of the auxiliary matrix A e
The values of the terms (A e i,j ) are given by
top we plot the function M (n, θ 0 ) for µ = 0.1, C = 6 and the values of n = 2 (continuous line) and n = 4 (discontinuous line).We remark its sinusoidal behaviour in accordance with equation (49).Consistently with Theorem 2, the curve M (n, θ 0 ) intersects four times M (n, θ 0 ) = 0.The four specific values of θ 0 give rise to four n-EC orbits belonging to families α n , β n , γ n and δ n .The corresponding EC orbits are shown in the bottom figure in usual synodical coordinates (x, y).

Figure 5 :
Figure 5: µ = 0.1, Top.We plot the angular momentum M (2, θ 0 ).Notice the zoom area where the appearance of two new bifurcating orbits (in green), besides the family α 2 is observed when decreasing C. Bottom.Left, middle and right.Four 2-EC orbits (the colour code corresponds to the top figure) for C = 3.76 (in blue), C bif = 3.72442505 (the bifurcating value, in red), C = 3.69 (in green).Darker colour: those 2-EC orbits belonging to family α 2 .In the right plot, also the two new bifurcated 2-EC orbits are plotted.

Figure 6 :
Figure 6: µ = 0.1, n = 3. Top.We plot the angular momentum M (n, θ 0 ).Notice the zoom area where the appearance of four new bifurcating orbits (in green), besides the family α 3 is observed when decreasing C. Bottom.Left, middle and right.Four 3-EC orbits (the colour code corresponds to the top figure) for C = 3.9 (in blue), C bif = 3.80644009 (the bifurcating value, in red), C = 3.7 (in green).Darker colour: those 3-EC orbits belonging to family α 3 .In the middle plot, also the two tangent new bifurcated 3-EC orbits are plotted.In the right plot, also the four new bifurcated 3-EC orbits are plotted.
]) we have computed the function M (n, θ 0 ) for n = 1, ..., 8.In Figures7 and 8we plot the obtained results for µ = 0.1 and µ = 0.8, what we call bifurcation diagrams.For n fixed, we plot the diagram (θ 0 , C) and the colour standing for the value of M (n, θ 0 ).The drastic change of colour (from yellow to green) describes the change of sign of M (n, θ 0 ) and therefore the existence of an n-EC orbit.So for any C fixed, we clearly see the number of n-EC orbits.Some comments about Figure7must be made: (i) for big values of C, the bifurcation diagrams show clearly four n-EC orbits for any value of n, in accordance with Theorem 2. See any plot in the figure.

Figure 11 :
Figure 11: Value of the angular momentum of the ejection orbits at the n intersection with Σ m for K ∈ [K L1 , 8] and n = 3, ..., 10.In black the values corresponding to n-EC orbits.

Figure 12 :
Figure 12: Top.Initial conditions for the 5-EC orbits corresponding to the families α n (yellow), β n (green), γ n (blue), δ n (red) and the new families of orbits (purple) as function of K. Bottom.The trajectories of the orbits (in correspondence with the previous color) that exist for the values of K denoted previously.The values of K correspond to the value of the bifurcation K ≈ 5.02714993 (left), a value where we have eight 5-EC orbits K = 4.86 (middle) and the value of collapse K ≈ 4.72835275.