Study of the ejection/collision orbits in the spatial RTBP using the McGehee regularization

. In this paper we analyse the McGehee regularization of the collision in the spatial Restricted three-body problem (3D RTBP). As a particular application, we study the ejection (collision) orbits. The parameterization method is applied up to high order to obtain suitable initial conditions of ejection (collision) orbits. Moreover, assuming ejection orbits, diﬀerent methods are discussed to detect which of them are also collision orbits. Finally we explore the so called ejection-collision (EC) orbits, that is, orbits where the particle ejects from one primary, reaches a maximum in the distance with respect to the same primary, and ends at collision with that primary. Some results concerning the existence of spatial EC orbits are described.


Introduction
As is well known the 3D RTBP consists in describing the motion of a particle subject to the gravitational attraction of two point massive bodies P 1 and P 2 , called primaries, assuming that the primaries describe circular orbits around their common center of mass.In a rotating system of coordinates and using suitable units of time, mass and length, the system of ODE to describe the motion of the particle is an autonomous nonlinear system of 6 ODE of first order, which has a first integral called Jacobi constant.
Many papers have been devoted to study the RTBP from a theoretical point of view.But also the RTBP has been used as a preliminary model to explain some real missions, for example Genesis and Artemis, or horseshoe motion of near Earth asteroids (see [1]), or the coorbital motion of Janus and Epimetheus (see [10]), or to construct new trajectories with low energy transfers (see [6]), to mention a few applications.Actually the planar RTBP (where the motion of the particle is restricted to the same plane of motion where the primaries are moving), has been more studied by far.But still there are some aspects that remain far from being well understood in both the planar and spatial cases.In particular, this paper is focussed on the EC orbits in the spatial RTBP.At present there are some results about such orbits in the planar RTPB (see [11,13,14,16] and references therein).The EC orbits in the spatial RTBP have been, however, hardly ever analysed.
A first important comment concerning the system of ODE for the RTBP (both planar and spatial cases) is that such a system has only two singularities that correspond to a collision with either one or the other primary.Since a main goal of this paper is to study EC orbits, that start/end at collision, it is clear that some kind of regularization must be applied.
In this framework, and focussing on the planar RTBP, there are several choices, from local ones (removing one of the two singularities -a collision with one primary-whereas the other one -a collision with the other one-remains) to global ones (regularizing both singularities).The McGehee's local choice is the easiest one from the point of view of the physical meaning of the regularizing variables considered (which are essentially polar coordinates); moreover, the system of ODE has a simple expression and one can analyse the collision manifold.However, ejection or collision orbits, in the McGehee variables (plus a scaling in time), become asymptotic solutions to equilibrium points and this is a drawback when doing numerical simulations.By contrast, Levi-Civita local regularization results more suitable from the numerical point of view, since the pass through collision is a regular point, but the expression of the resulting ODE is more intricate and a double covering of the phase space intrinsically appears.See [15] where a comparison between both choices is carried out and some other regularizations are also mentioned.
Once we move from the planar to the spatial RTBP, the first natural motivation is to consider the generalization of McGehee regularization (called from now on generalized McGehee regularization) into the 3D case.It is still true that the physical meaning is easy because we have to consider simply spherical coordinates (instead of polar ones in the planar case) plus a scaling in time.However following the same steps done to apply McGehee regularization to the planar RTBP, now in the spatial case, one obtains a system of ODE which is singular along the z axis.We remark that the z axis is not singular in the real, physical, problem; the singular axis appears simply by the choice of the spherical coordinates.
Notice that the generalized McGehee regularization is the first natural generalization when moving from the planar to the spatial RTBP.The next step is to analyse the generalization of Levi-Civita regularization which becomes Kustaanheimo-Stiefel's [18].In this case, however, besides the loss of physical meaning of the new variables, the regularization must deal with 4D spaces instead of the 3D natural space.
So a first goal of this paper is to focus on the generalized McGehee regularization for the spatial RTBP.Since the obtained system is not a suitable one to use numerically when we consider trajectories with high inclinations (passing close to the z axis requires big ranges of time), some strategy must be applied.Our approach has been to take two different local charts where the numerical integrations work perfectly well.
As a particular application of the system of regularized ODE obtained with the generalized McGehee regularization, we focus on the second goal of this paper which is the study of ejection (collision) orbits, and also ejection-collision (EC) orbits.Since ejecting from (colliding at) a primary means leaving (reaching) an equilibrium point along an asymptotic trajectory -that is the unstable (stable) invariant manifold of the equilibrium point-, the parameterization method has been applied in order to obtain suitable initial conditions of such trajectory.We notice that by means of the parameterization method, typically expansions in some parameter s close to the origin are obtained to have a parameterization of an invariant manifold.However, due to the expression of McGehee system of ODE, the parameterization is an expansion in terms of s 1/2 .Once we know how to eject from (collide at) a primary, we focus specifically on the EC orbits.We carry out numerical simulations in order to visualize the existence of EC orbits in the spatial RTBP.At this point it is important to mention the unique reference, [9], we know about the EC orbits in the spatial RTBP.In that paper, the authors prove the existence of two circles of EC orbits, for µ > 0 fixed and a value of the Jacobi constant C ≥ C L1 fixed (being C L1 the value of the Jacobi constant at the collinear equilibrium point L 1 ).We also recall that the condition C ≥ C L1 guarantees that there is no possible transit from moving in a neighborhood of one primary to a neighborhood around the other one.Concerning the results from our simulations, we will show how their statement about the existence of two circles of EC orbis is not true, but only six EC orbits can be guaranteed for µ > 0 and C ≥ C L1 fixed: four of them being planar and already known (see [11,14,17]) and two new spatial EC orbits.
In short, the main contribution of this paper is twofold: (i) on the one hand, we consider the McGehee regularization in the spatial RTBP, we explain our approach to deal with highly inclined trajectories and we apply the parameterization method to efficiently find the initial conditions of ejection (collision orbits.(ii) On the other hand, we analyse the existence of EC orbits in the spatial case.
The paper is organized as follows.In Section 2 we recall some well known properties of the spatial RTBP that will be used along the paper.In Section 3 we apply McGehee regularization and introduce the new scaling in time.In Section 4 we explain the strategy to deal with orbits with high inclination.The strategy is based on taking two (out of three) different local charts.Section 5 is devoted to describe shortly the flow on the collision manifold.In Section 6 we explain the parameterization method to obtain initial conditions of ejection (collision) orbits.Section 7 is focussed on the numerical strategy to compute EC orbits.Finally, we provide the results in Section 8.
We remark that all the computations have been done using double precision and the numerical integration of the systems of ODE has been done using an own implemented Runge-Kutta (7)8 integrator with an adaptive step size control described in [5] and a Taylor method implemented on a robust, fast and accurate software package in [8].

The RTBP
The circular, restricted three-body problem (RTBP) describes the motion of a particle of infinitesimal mass, moving under the gravitational influence of two massive bodies, called primaries, that describe circular orbits around their common center of mass on a given plane.We will consider the spatial (or 3D) RTBP, in which the motion of the particle is allowed on the whole space of configuration.Taking a coordinate system that rotates with the primaries, with origin located at their center of mass, and suitable units, we can assume that the primaries have masses 1 − µ and µ, µ ∈ (0, 1), their positions are fixed at (µ, 0, 0) and (µ − 1, 0, 0), respectively, and the period of their motion is 2π.With these assumptions, the equations of motion for the particle in this rotating (also called synodical system) are given by where and ˙= d/dt.
It is well known that this system of ODE has the following properties (see [19] for details) which will be used along the paper: 1.It has a first integral, the so-called Jacobi integral, defined by 2. The equations of motion satisfy the symmetries This means that, given any solution on the configuration space (x(t), y(t), z(t)) of system (1), there exist symmetrical ones with respect to the (x, y) plane and with respect to the (x, z) plane.These symmetries will be useful later on for our purposes.
3. There exist 5 equilibrium points (with ( ẋ, ẏ, ż) = (0, 0, 0)): the collinear ones, L i , i = 1, 2, 3 on the x axis, and the triangular ones L i , i = 4, 5 located on the (x, y) plane at the vertices of an equilateral triangle with the primaries.We will assume that , that is, L 1 is between the primaries, L 2 is on the left hand side of the small one and L 3 on the right hand side of the big one.We will denote by C Li (µ) the value of the Jacobi constant at L i for a given µ.
4. The equations of motion can be written as a Hamiltonian system in position (x, y, z) and momenta (p x , p y , p z ) variables, with p x = x − y, p y = y + x and p z = z defined by the Hamiltonian function with r 1 = (x − µ) 2 + y 2 + z 2 and r 2 = (x − µ + 1) 2 + y 2 + z 2 , and the relation between C and H is given by We denote by H Li (µ), the associated value of the Hamiltonian at L i for a given µ. 5. Fixed a value of the Jacobi constant C (or the Hamiltonian H), the motion is allowed to take place in the Hill's region defined by In this paper we will restrict the range of values of C to C ≥ C L1 (µ) (H ≤ H L1 (µ)).As we can see in Figure 1, for such a range of values of C the bounded motion can take place only in a region around each primary (the red and purple regions on the plot) and there is no possibility of trajectories going from one to the other primary.Of course, we might consider ejection/collision orbits for smaller values of C, but in such a case other invariant objects appear (periodic orbits, homoclinic and heteroclinic ones, 2D tori, etc) which interact with the ejection orbits and the dynamics is very rich and much more intricate.See a (partial) description of this dynamics for the planar RTBP in [16].Along the paper, we will use specific values of H, that can be translated to values of C through relation (6).

McGehee regularization in the 3D case
In order to study orbits that eject from or collide with a primary, we need to apply some regularization strategy, since the system (1) and also the Hamiltonian ( 5) are singular at the collision (that is r 1 = 0 or r 2 = 0).In this paper, we want to analyse the McGehee regularization, which is local in the sense that, after applying some changes of variables and time, we obtain a new system of ODE which is regular at r 1 = 0 but singular at r 2 = 0 (or regular at r 2 = 0 and singular at r 1 = 0).In the sequel, we describe the different steps to regularize the collision r 1 = 0.
The first step is to move the primary P 1 to the origin, by means of a translation: x = x − µ, ȳ = y, z = z, px = p x = ẋ − y = p x , py = p y = ẏ + x = ẏ + x + µ, pz = p z .We obtain a new Hamiltonian, we abuse the notation dropping the bar and denoting by (x, y, z, p x , p y , p z ) the new variables.The Hamiltonian of the 3D RTBP is: where The next step is to consider the symplectic change to spherical variables given by with the longitude θ ∈ [0, 2π) (on the (x, y) plane) and the inclination (or latitude ].We will call this change spherical variables with respect to the z axis. Let us recall that the symplectic change of variables (q, p) to (Q, P ) in R 2n , when we know q = f (Q) comes from: In our case, we have p = (p x , p y , p z ) and P = (p r , p θ , p ϕ ).So, we have and the Hamiltonian (7) becomes: with r 2 = r 2 + 2r cos θ cos ϕ + 1.
The associated system of ODE is given by From (10) and (11) we have: As it is done in the 2D case, we introduce the variables: and a rescaling in time through the relation With these changes system (10) becomes At this point, we want to emphasize that if we follow this same strategy as in the 2D case, we still have a problem: the resulting system of ODE is not defined for ϕ = ±π/2, i.e. all the points on the z axis are singular.This situation did not happen in the 2D case and it is precisely the reason why the regularization of Kustaanheimo-Stiefel requires an extra dimension (see [18] for details).Actually the singularity introduced along the z-axis does not correspond to a physical singularity, but is caused by the spherical coordinates and can be removed considering the change of time dt/dτ = r 3/2 cos ϕ instead of the change (13).With this new change of time only the singularity r 2 = 0 persists; however the z axis (cos ϕ = 0) -on the configuration space-is invariant.For this reason, both the system with this new change of time or system (14) are still not convenient systems to work with orbits that pass through or near the z axis and other strategies are preferable as we will see in the next section.
We also remark that the first integral (the Jacobi constant or the Hamiltonian itself) is defined neither on r = 0 nor on r 2 = 0, However if we multiply the relation H = H by r we obtain which is valid on each energy level H = H constant.

Local charts
From the above discussion it is clear that one drawback appears: the system ( 14) is not a suitable one to deal with numerically if we want to consider trajectories with high inclinations (that is passing close to the z axis), since the error accumulated may become very large.
Our approach to deal with this inconvenience consists in taking two (out of three) different local charts in such a way that, for each chart we only consider inclinations in the interval [−π/4, π/4] -and the numerical integrations work perfectly fine-, and the union of the local charts recover the whole space.
More precisely, we consider a first chart taking into account system (14), that is using the McGehee variables from the spherical variables with respect to the z axis.We avoid high inclinations just   We emphasize that any pair of local charts suffices to cover the whole sphere S 2 .In spite of this we consider the other two changes of spherical variables and proceed to obtain the corresponding systems of regularised ODE.We simply provide the resulting expressions.

Spherical variables with respect to the x axis
We consider the symplectic change of variables to spherical variables with respect to the x axis, that is A similar procedure, with = d/dτ where now dt/dτ = r 3/2 , leads to the system where r 2 = r 2 + 2r sin ϕ + 1.

Spherical variables with respect to the y axis
We consider the symplectic change of variables to spherical variables with respect to the y axis, that is,      x = r sin θ cos φ, y = r sin φ, z = r cos θ cos φ.
Remark It is clear that using a combination involving two of the systems ( 14), ( 17) and ( 18), we can integrate any orbit, except when it passes very close to the second primary.In such a case we can proceed in a similar way with the same kind of regularization around the second primary.

The collision manifold
The collision manifold, called from now on Λ, is the manifold invariant by the flow of system (14) (or (17) or ( 18)) when we take r = 0. We will analyse the collision manifold in variables (r, θ, ϕ, v, u θ , u ϕ ).
A similar study can be done when taking the other set of variables introduced for the different local charts.
A first observation is that the energy relation ( 16) on r = 0 becomes u 2 ϕ + u 2 θ + v 2 = 2(1 − µ) for arbitrary θ and ϕ ∈ (−π/2, π/2).Therefore, considering and extra local chart we have that the manifold Λ is a four-dimensional manifold diffeomorphic to S 2 × S 2 and is defined on the boundary of each fixed level of the Jacobi integral (or the Hamiltonian).In particular, the system of ODE ( 14 which describes the internal dynamics on the manifold.
On Λ there exist two spheres S 2 ± of equilibrium points defined by In order to determine their stability as equilibrium points in R 6 , we proceed to the linearization of system ( 14) at any point belonging to S 2 + with ϕ ∈ (−π/2, π/2) and we obtain the Jacobian matrix: Similarly we obtain the linearization of the system on each point of The corresponding eigenvectors are: Let us denote by W u,s (P ) the unstable and stable invariant manifolds of an equilibrium point P , and W u,s (S 2 ± ) the unstable and stable invariant manifolds of the whole sphere S 2 + or S 2 − .Concerning the dimensions of these invariant manifolds and the internal dynamics on Λ, Llibre and Alfaro in [9] proved the following result.
Proposition (i) Fixed an energy level we have (ii) For the flow on the collision manifold Λ, the unstable invariant manifold associated with the equilibrium point , coincides with the stable invariant manifold of the equilibrium point This can be intuitively understood as follows: Once an energy level is fixed, each point P + of S 2 + imposes the ejection direction from the first primary (similarly each point Q − of S 2 − imposes the collision direction).Taking into account the time, fixed an energy level we have that W u (S 2 + ) is diffeomorphic to S 2 × R and therefore its dimension is 3.
The other relation is obtained directly from the internal dynamics of the collision (19).The collision manifold lives in every energy level and, as we have said previously, is diffeomorphic to S 2 × S 2 .In particular (see [9] for details) we have that this collision manifold is foliated by heteroclinic connections from points of S 2 − to points of S 2 + .
At this point we distinguish between 3 types of orbits that start/end (asymptotically in time) in the collision manifold, r = 0, but that along the trajectory r > 0 and their study will be the purpose of the remaining sections: (i) ejection, (ii) collision and (iii) ejection-collision orbits.More specifically: Definition 1.The ejection orbits manifold, W u (S 2 + ) is the set of ejection orbits -those for which the particle ejects from a primary-, that is, the set of orbits on the unstable manifold W u (P + ), for any P + = (0, θ, ϕ, v 0 , 0, 0) ∈ S 2 + .So each ejection orbit may be regarded as an orbit such that r > 0 for all finite time τ and asymptotically tends to an equilibrium point P + ∈ S 2 + as τ → −∞.
Definition 2. The collision orbits manifold, W s (S 2 − ) is the set of collision orbits -those for which the particle arrives at collision with a primary-, that is, the set of orbits on the stable manifold W s (Q − ), for any Q − = (0, θ, ϕ, −v 0 , 0, 0) ∈ S 2 − .So each collision orbit may be regarded as an orbit such that r > 0 for all finite time τ and asymptotically tends to an equilibrium point Q − ∈ S − as τ → +∞.Definition 3. The set of ejection-collision orbits, EC orbits, -those for which the particle ejects from/arrives at collision with the same primary-is the set of orbits obtained from the intersection − ).So they may be regarded as heteroclinic orbits between In order to somewhat classify the EC orbits, another definition existing in the literature is the so called n-ejection-collision orbit: Definition 4. 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.
These kind of orbits have been studied in depth only for the planar case (see for example [11,13,14,15,16,17]).
We show in Figure 4 examples of 1-EC orbits in the planar RTBP for µ = 0.1, 0.5.Actually, in this paper we focus on 1-EC orbits for the spatial RTBP.

Ejection/collision in the McGehee regularized system
In order to compute numerically the set of ejection (collision) orbits, we need to compute the unstable (stable) manifold of any equilibrium point belonging to S 2 + (S 2 − ).To do so, fixed a value of H = H, in order to consider initial conditions for an ejection orbit, we make two comments: on the one hand, for each equilibrium point P + ∈ S 2 + , the 2D W u (P + ) is tangent to the plane passing through P + generated by the eigenvectors v 1 and v 2 of (21), i.e. vectors like v = (a, 0, 0, b, 0, 0), with a, b ∈ R. On the other hand, we recall that the energy level set H = H is defined implicitly by (16), and the normal vector to this energy level set at point P + = (0, θ 0 , ϕ 0 , v 0 , 0, 0) is n = (−H − 3 2 µ, 0, 0, v 0 , 0, 0).So the vectors v (generating the tangent plane W u (P + )) must be perpendicular to n, i.e. they must satisfy In this way, the normalized vector tangent to W u (P ) for H = H is given by As we want to consider the linear approximation to this manifold we will take the initial condition of an ejection orbit associated with the point P = (0, θ 0 , ϕ 0 , v 0 , 0, 0) as with s > 0 a small quantity (typically 10 −6 to 10 −8 ).Varying θ 0 and ϕ 0 and considering another local chart we can cover the whole set of initial conditions of ejection orbits.This procedure has the main advantage of being very simple but it has a clear disadvantage: for obtaining a small error we must use a sufficiently small s and therefore we will need more integration time, which will cause numerical error to accumulate.In order to avoid this problem we will use approximations of high order that can be computed via the parameterization method (see [2,3,4]).

The parameterization method
The parameterization method allows us to obtain a high-order approximation of invariant manifolds associated with a dynamical system z = F (z).In our case we want to obtain the invariant manifold associated with the equilibrium point P + , tangent to the vector w 1 .
To simplify the notation, we will define z = (r, θ, ϕ, v, u θ , u ϕ ) and we rewrite the system ( 14) as z = F (z).To find a parameterization z = W (s) of the invariant manifold W u (P + ) associated with the energy level H = H (this is tangent to w 1 in P + ) such that W (0) = P + we have to solve the invariance equation where the internal dynamics of the manifold is given by the equation s = f (s) with f (0) = 0.
In order to solve the equation ( 26) the typical strategy is to find power series in terms of s of W and f (see [7]).However, due to the expression of the system (14), our case does not admit a series expansion in s since we want to expand at points with r = 0 and in the system we have terms with r 3/2 .To deal with it we must consider an expansion in powers of s 1/2 .
To simplify the problem to be solved, we consider the system where S is the matrix that has column vectors the eigenvectors of DF (P + ) given by w 1 and , 0, 0, −v 0 , 0, 0, , where the expression of the vectors v i , i = 3, ..., 6 is given in (21).In this way the equilibrium point of ( 27) is 0 and and the equation that we have to solve becomes with W (s) = P + + SW(s) that implies W(0) = 0. Thus we have to find W(s) and f (s) as an expansion in terms of with ω 1 = (1, 0, 0, 0), and we will introduce subscript <k/2 or ≤k/2 in order to denote the truncation of the expansion until the order (k − 1)/2 and k/2 respectively.For example At this point the goal is to compute the terms ω k/2 and f k/2 for k > 1 once we know W <k/2 (s) and f <k/2 (s) (i.e.knowing W 1 , ... , W (k−1)/2 and f 1 , ... , f (k−1)/2 ).
We know the expression on the left of the invariance equation (29) up to order (k − 1)/2, i.e.

G(W
. Thus, we will first compute the homogeneous terms of degree k of the composition G(W <k/2 (s)) and of the matrix product DW <k/2 (s)f <k/2 (s) which we will denote by: In this way, the terms of order k/2 of the invariance equation ( 29) lead us to the cohomological equation of order k/2 for W k/2 and f k/2 where Therefore, considering a normal form parameterization style we have where η j k/2 denotes the j = 1, ..., The first terms of the manifold W (s) are given by where h = H + 3µ/2 and B = h 2 + v 2 0 .
It is important to note that s ≥ 0, since r ≥ 0 and that we have the 2D parameterization of the manifold W in terms of (H, s).This procedure to compute an approximation of the manifold can be done up to the desired order, but in practice we have usually worked with the approximation of order 5-10 in s and with values of s in the range 10 −3 to 10 −2 .
An important remark at this point is that, we can use a similar procedure in order to compute the approximation of the manifold using the local charts described in the previous Section, the whole set of initial conditions of ejection (collision) orbits (that is for all value of the longitude θ ∈ [0, 2π) and latitude ϕ ∈ [−π/2, π/2]) can be considered.So with this procedure we obtain the whole set of ejection orbits, W u (S 2 + ) and similarly, we can compute the whole set of collision orbits, W s (S 2 − ).
In summary, we have described how to compute the whole set of initial conditions of ejection orbits and of collision ones.The next natural step is to analyse how to detect the existence of collision orbits, and, more particularly, EC orbits.This is the purpose of the next Section.
7 Detection of candidates to collision orbits.EC orbits.
In order to detect a collision we can use different strategies (see [17]).In particular, two suitable options are based on either checking that a point belongs to an orbit of the collision orbits manifold or on the computation of the angular momentum at the minimum distance from the primary.We describe both approaches and discuss the pros and cons when applied to the planar and spatial RTBP and the computation of EC orbits as well.

Orbits that belong to the collision orbits manifold
Let us discuss a first strategy to guarantee the existence of collision orbits based on the intersection of the collision orbits manifold with a Poincaré section Σ.The idea is simply to prove that there exists (at least) one point in Σ that belongs to the collision orbits manifold.
Let us compare this strategy for the planar and spatial RTBP.Given a value of µ ∈ (0, 1), a main difference between both problems is the following: in the planar case, we have four dimensions, and a 3-dimensional space when a value of the Jacobi constant is fixed, and a 2-dimensional one when a Poincaré section is taken.The collision orbits manifold is 2-dimensional, so once we integrate it backwards in time and look at the intersection between the manifold and a Poincaré section, we will have a curve C s which is 1-dimensional.On the other hand, assume we have a (piece of) curve α, corresponding to the intersection between a set of orbits -that may be ejection orbits or not-and the Poincaré section.Any point belonging to α ∩ C s integrated forward in time is a collision orbit.Therefore playing with a transversal segment with 2 end suitable points (say outside and inside the curve C s respectively), we can apply Bolzano's theorem to detect the existence of an intermediate point on the curve that is precisely a collision orbit.
However in the spatial RTBP, we start with 6 dimensions and reduce two dimensions (fixing a value of the Jacobi constant and a Poincaré section).The collision orbits manifold is 3-dimensional and the intersection with the Poincaré section gives rise to a surface, that is a 2-dimensional manifold in a 4-dimensional space.Since we do not have the inside and outside notion, we need a suitable transversal 2-dimensional domain and it is not that easy to guarantee the existence of collision orbits.However this approach is not straightforward in the spatial RTBP.As previously mentioned now the ejection orbits manifold (W u (S 2 + )) and the collision orbits one (W s (S 2 − )) are 3-dimensional and the intersection with Σ M become two surfaces (2-dimensional manifolds) in a 4-dimensional space.It is not so clear how to detect the intersection points between the two surfaces and we will show precisely this difficulty with some simulations in Section 8.

Angular momentum at the minimum distance from the primary
The second method that we consider to detect collision orbits is based on the angular momentum, or quantities related to it.Assume a value of µ fixed and C given.Take a small neighborhood centered at the origin and of radius > 0 small.Then the method to detect collision orbits is based on this statement: the conditions where m is the angular momentum vector, are not compatible.Indeed, from the two first conditions the particle reaches a point Q = (x, y, z) or (r, θ, ϕ) which is a minimum in distance with respect to the origin.If > 0 is small enough, the motion of the particle will be essentially the motion of a two body problem and the velocity will be a non zero vector (the modulus of the velocity can be obtained from 2Ω(x, y, z) − C = ||v|| 2 and the point Q is far from the zero velocity surface).At Q, with r = 0, the position velocity and vector velocity will be orthogonal vectors, and therefore it is impossible for the angular momentum vector to be equal to zero.Therefore, for > 0 small enough, the requirement of the first three conditions implies that r = 0, that is a collision between the particle and the primary.This is precisely the numerical strategy we have implemented to compute collision orbits.
We now discuss how this strategy works very well in the planar RTBP but not in the spatial case.

Detection of EC orbits in the planar RTBP
Let us shortly explain this strategy to detect EC orbits in the planar RTBP.It is based on the change of sign in the angular vector at a minimum distance with respect to the primary the particle ejected from.This means that one single condition must be satisfied.(This will not be the case when the same strategy is applied to the spatial RTBP.)Consider now the planar RTBP.In order to detect an EC orbit, we start at ejection orbits, and for each one we compute the angular momentum vector at the minimum distance with respect to the origin.An EC orbit will be that one such that the angular momentum at that minimum will be zero.Of course, we cannot reach exactly an EC orbit, because it is an heteroclinic connection between two equilibrium points and the integration time is infinite.Since the angular momentum vector is m = (0, 0, m z ), our numerical strategy will consist of looking for a change of sign in m z .This will be a transversal condition very suitable for numerical simulations.We proceed as follows: we take the whole set of ejection orbits parametrized by θ 0 (see ( 33)) (and using the parameterization method described in Section 6 for the planar RTBP) and integrate them up to the first path to minimum distance with the primary they ejected from, that is we compute the intersection of the set of ejection orbits up to the first crossing with the Poincaré section Σ m , assuming that this distance is small.At such intersection we compute m z and we plot the curve (θ 0 , m z ).The initial condition (corresponding to the value of θ 0 ) of an EC orbit will be detected at each crossing of the angular momentum curve with m z = 0. Let us show this argument geometrically in Figure 6, middle and right.On the middle plot, the purple and yellow ejection orbits are plotted in continuous line from ejection to the first minimum in the distance r and in discontinuous line for a small range of time after the minimum distance.We clearly see that the signs of the corresponding angular momentum vectors (at the path of minimum distance r) for the purple and yellow ejection orbits are opposite; therefore by continuity, there must exist a suitable value of θ 0 such that at the first minimum distance the angular momentum is equal to zero, that is, the ejection orbit has a collision with the primary, and it becomes an EC orbit (the black trajectory in the middle plot).We also show this argument in Figure 6 right.The green and orange lines correspond to the heteroclinic connections between P + and P − in Figure 6 left.The purple and yellow trajectories move at different sides close to P − , the purple one passes through a point at Σ 0 with v = 0 and m z > 0 (or u θ > 0) and the blue one through a point at Σ 0 with v = 0 and m z < 0 (or u θ < 0).So by continuity there is an intermediate ejection orbit that will tend to a point Q ∈ S 1 − , that is m z will be zero (u θ = 0) and this is an EC orbit.See the black piece of curve tending to P − on the right plot.
So, in summary one single transversal condition must be satisfied: the curve (θ 0 , m z ) must cross the axis m z = 0 assuming that the ejection orbits close to the EC one have a small value of r at the first crossing with Σ m .It is clear that eventually an ejection orbit might have a value of m z = 0 at the first crossing with Σ m but a big value of r thus not being an EC orbit.Actually, numerically, we check that the ejection orbits really pass close to the primary at the first crossing with Σ m .The interested reader is referred to [11,13,14,16,17] for massive numerical simulations on n-EC orbits, for n ≥ 1, in the planar RTBP.

Detection of candidates to EC orbits in the spatial RTBP
Now let us discuss the strategy based on the angular momentum in the spatial RTBP.We recall that on the collision manifold (r = 0) there is a whole sphere S 2 + of equilibrium points and a whole sphere S 2 − of equilibrium points.Moreover there is a sphere connecting asymptotically the point − with the point P + = (0, θ 0 , ϕ 0 , v 0 , 0, 0) ∈ S 2 + , that is, we have an sphere of heteroclinic connections between P − and P + .See the grey curves in Figure 7 left.We recall that this sphere of grey curves simply becomes the union of two heteroclinic connections in the planar RTBP (see Figure 6).First let us discuss how to detect a collision orbit.Now the geometric picture is the one described in Figure 7. Roughly speaking, in the spatial problem we need a circumference of trajectories (instead of two orbits in the planar case) such that cover the whole grey sphere, forward in time.In such a case there is an orbit that tends to P − and thus, is a collision orbit.See Figure 7 left and compare with Figure 6 right.We also compare the plot in Figure 6 middle for the planar case with the ones in Figure 7 middle and right.We take now a cone of ejection orbits in the (x, y, z) configuration space (plus the interior ejection orbit) and follow the evolution of the orbits up to the first passage to minimum distance (with respect P 1 located at the origin).The whole evolution from ejection to close to collision is globally shown in the middle plot.On the right plot we see a zoom region close to the origin: the cone of ejection orbits (plus the black one inside that will be an EC orbit) and the kind of umbrella set close to collision (plus the collision orbit in black inside).The single black orbit is an EC orbit.Now let us focus on the numerical strategy to detect EC orbits based on the angular momentum vector which now is m = (m x , m y , m z ).We take the whole set of initial conditions of ejection orbits whose initial conditions are parametrized by θ 0 and ϕ 0 .In principle ϕ 0 ∈ [−π/2, π/2], but as described above, we will consider different local charts.Our goal is to obtain a value (or several values) of (θ 0 , ϕ 0 ) such that the angular momentum at the crossing with Σ m , is equal to zero.Such (θ 0 , ϕ 0 ) will provide the initial condition of an EC orbit.
So, for each (θ 0 , ϕ 0 ), we follow the ejection orbit up to Σ m .Rather than in m, we are interested in its direction defined by two angles (α, δ) through The actual numerical computation of the angles α and δ follows from two steps.First the angular momentum m = (m x , m y , m z ) for each ejection orbit at the crossing with Σ m is computed taking spherical variables with respect to the corresponding axis considered and in terms of the McGehee variables and time.For example, taking spherical variables with respect to the z axis, we obtain Similarly, we obtain the corresponding expressions if the other sets of spherical variables are used.Second, we compute the associated angles α and δ.
• (ii) Taking into account α, we look for a sudden change of complete direction of the vector m in the (x, y) plane, that is, α → α + π.We obtain a curve or set of curves labelled by A.
• (iii) Taking into account δ, we look for a sudden change of complete direction of the vector m in the (z) direction, that is, δ → −δ.We obtain a curve or set of curves labelled by D.
• (iv) The values of (θ 0 , ϕ 0 ) such that the curves A and D intersect correspond to initial conditions of EC orbits.

Numerical results on 1-EC orbits
We describe now the numerical results obtained for the 1-EC orbits, using the two different methods discussed in the previous Section.
Applying the first method, we consider the Poincaré section Σ M (maximum distance with respect to the origin).We integrate the set of ejection orbits (forwards in time) up to the first crossing with Σ M and similarly we integrate (backwards in time) the set of collision orbits up to the first crossing with Σ M , obtaining the associated points (r, θ, ϕ, v, u θ , u ϕ ) u and (r, θ, ϕ, v, u θ , u ϕ ) s , respectively.We must look at the intersection points.Notice that since v = 0 and a value of the Jacobi constant is fixed, we just need to check that four out of six variables are equal when comparing the ones coming from the ejection and the collision.So we consider two sets of three variables (θ, ϕ, u θ ) u,s and (θ, ϕ, u ϕ ) u,s .They are shown in Figure 8 in red and blue colors respectively.Although it is true that in the left figure we see that the surfaces intersect transversally in four curves, we notice in the right one the intersection is along two curves that are (seem to be) tangent.Actually, even in the case that the intersection were transversal, we should need an additional check to guarantee that the coincident curves from the left and right plots take place exactly for the same values of θ and ϕ.Now let us apply the second method.We take the Poincaré section Σ m , and we follow the four steps (i),...,(iv) explained above.For a fixed µ ∈ (0, 1) and C ≥ C L1 given, we integrate just the whole set of ejection orbits, as before, up to the first crossing with Σ m .Now we are interested in the values of α and δ of the angular momentum vector m, obtained in the corresponding set of variables.We will call them the local parametrization of α and δ, respectively.So we have two local parametrizations for each local chart.
In Figure 9 we plot the resulting local parametrizations for µ = 0.1 and C = 4.We plot the value of α(θ 0 , ϕ 0 ) on the top plots and the value of δ(θ 0 , ϕ 0 ) on the bottom ones.For clearer visualization, we also plot the values of α and δ on the sphere S 2 with the usual variables (x, y, z) and the usual longitude and latitude θ and ϕ.See Figure 10.
In Figure 11 we plot the whole composition plot α(θ 0 , ϕ 0 ) and δ(θ 0 , ϕ 0 ).This corresponds to step (i).Now we proceed with steps (ii) and (iii).We plot the resulting curves A and D. The changes α → α+π, result in two curves A 1 and A 2 , that is A = A 1 ∪ A 2 ; and the changes of sign in δ result in two curves D 1 and D 2 , that is D = D 1 ∪ D 2 .They are shown in blue and red respectively in Figure 12.

Finally step (iv) requires the intersection points belonging to
Concerning the intersection of both sets of curves, we distinguish: • Four transversal intersection points with ϕ 0 = 0. See the two blue points and two red ones in Figure 13 left.The associated values of θ 0 provide the initial conditions of four 1-EC orbits which are planar.These four planar 1-EC orbits were expected in accordance with the analytical results proved in [11,14,17].See also Figure 4.The two blue points in Figure 13 left provide the initial conditions of the two 1-EC orbits which are symmetrical with respect to the x axis    also in blue in Figure 4.The two red points in Figure 13 provide the initial conditions of the two 1-EC orbits which are non symmetrical, but symmetrical one of the other one with respect to the x axis, also in red in Figure 4. • Two transversal intersection points with ϕ 0 = 0.The associated values of (θ 0 , ϕ 0 ) provide the initial conditions of two 1-EC orbits which are non planar.See the green dots in Figure 13 left.Each 1-EC orbit is itself symmetrical with respect to the plane xz, and symmetrical to the other orbit with respect to the plane xy.
• Apparently there appears to be a whole continuous family of intersection points that belong to A 2 ∩ D 2 (see Figure 12), plotted in grey discontinuous line in Figure 13 left.We cannot guarantee numerically that the curves A 2 and D 2 coincide so we can claim that there appears to be a continuous family of candidates to non planar 1-EC orbits.Nevertheless, we have done the following numerical computation: for each point (θ 0 , ϕ 0 ) we have followed the corresponding ejection orbit up to the first minimum distance from the origin (where the primary is located) obtaining the plot in Figure 13 right.We clearly see very small values of such distance (order of 10 −12 , but we might reach smaller values) precisely on the points belonging to the curves A 2 and D 2 , apart from the four points belonging to both curves that correspond to 1-EC orbits.We recall that given any 1-EC orbit, due to the symmetries (4), there exist three additional 1-EC orbits, so the hypothetical 1-EC orbits in the candidate family must also satisfy this property.
Summarizing, for µ > 0 and C ≥ C L1 fixed, our results are the following: • there are two planar 1-EC orbits which are symmetrical with respect to the x axis.Their initial conditions are provided by the values (θ 0 , ϕ 0 ), with ϕ = 0, in blue in Figure 13 for µ = 0.1 and C = 4.
• There are two planar 1-EC orbits which are non symmetrical with respect to the x axis, but symmetrical one of the other one with respect to this axis.Their initial conditions are provided by the values (θ 0 , ϕ 0 ), with ϕ = 0, in red in Figure 13.
• There are two spatial 1-EC orbits, symmetrical themselves with respect to the (x, z) plane, and symmetrical one of the other one with respect to the (x, y) plane.Their initial conditions are provided by the values (θ 0 , ϕ 0 ) in green in Figure 13.We have computed both the apparently coincident curves and the four points (that provide the initial conditions of 1-EC orbits) varying both C and µ.We plot the results obtained in Figure 14 for µ = 0.1 (top) and for µ = 0.5 (bottom), taking values of C = C L1 , 4, 6, 8, 10, 20, 50, 100.For µ fixed, we plot the curve and the four points corresponding to 1-EC orbits with the same colour.We observe that when C increases, the values of θ 0 tend to 0, π/2, π and 3π/2 for the planar 1-EC orbits.These results are in accordance with the known properties of planar 1-EC orbits (see [11,14,17]).Concerning the two spatial 1-EC orbits, we observe that as C increases, (θ 0 , ϕ 0 ) tends to (π, ±π/2).
From the numerical results obtained, we can claim that, for µ fixed and varying C ≥ C L1 , we obtain four families of planar 1-EC orbits (as already known) and two new families of spatial 1-EC orbits.Moreover, for µ fixed and C ≥ C L1 fixed, a continuous family of candidates to 1-EC orbits also does exist.
Notice that, as expected for µ fixed, as far as C increases, the maximum distance of the candidates to 1-EC orbits decreases, since the bounded Hill regions of possible motion around each primary (for C ≥ C L1 ) decrease.This effect is clearly shown in Figure 15, comparing the left plot for C = C L1 = 3.68695322987989 (left) and C = 6 (right) where the grey cloud represents the whole set of candidates to 1-EC orbits in the (x, y, z) configuration space, apart from the 1-EC orbits plotted in red, blue and green.Finally, we make a remark about the paper [9].In that paper the authors claim there exist, for small µ > 0 and fixed C > C L1 , two families of spatial, symmetric themselves (with respect to the xz-plane) 1-EC orbits, parametrized by the angle ϕ.For µ = 0 every ejection orbit ends in collision and an argument based on the symmetries of the problem is used to show that some of them survive when µ > 0. It can be easily seen that this is indeed the case for two planar (ϕ = 0) 1-EC orbits.The symmetry argument given in [9] to show the existence of a family of spatial 1-EC orbits (parametrized by ϕ) is, however, unclear.Furthermore, our numerical results contradict the existence of such a family or families: the planar symmetric themselves 1-EC orbits (the blue ones in Figure 15) are isolated and do not continue to a family (or families) of spatial 1-EC orbits.This non existent family (families) is precisely the one the authors there claim to exist, but our numerical results allow to conclude that they do not.As shown in Figure 15, the family of candidates to 1-EC orbits are continued from the red (non-symmetric) ones.

Conclusions
We have shown how the analysis of n-EC orbits, for n = 1 is more difficult in the spatial RTBP rather than in the planar case.We have considered just the case n = 1 to illustrate the difficulties.The natural regularizing McGehee coordinates plus the scaling in time, although natural and easier from the physical point of view, leads to a system of ODE that requires some additional treatment when doing the numerical simulations with high inclination orbits.This specific treatment is done using local charts.Moreover the analysis of existence of 1-EC orbits in the spatial RTBP is more tricky due to the dimensions involved.
Despite these complexities, the use of the parameterization method allows us to compute the initial conditions of the ejection (collision) orbits in an effcient way.The integration of these orbits from a point close enough to the ejection (or collision) is avoided with the parameterization method and therefore, the numerical cost of this strategy is reduced drastically.
It seems clear that in order to analyse the n-EC orbits, for n > 1, the KS regularization, although much less intuitive and that involves a system of eight (instead of six) ODE, is more convenient, since we have to deal with previous collisions and/or close encounters.This will be done in a future work.

Figure 1 :
Figure 1: µ = 0.2.Hill region of motion for C > C L1 (left), C = C L1 (right).The forbidden region of motion is limited between the red and purple surfaces and the grey one.
taking into account values of ϕ ∈ [−π/4, π/4].Similarly we introduce a second local chart using the McGehee variables from the spherical variables with respect to the x axis.We avoid high inclinations just taking into account inclination values in [−π/4, π/4].Finally we introduce a third local chart using the McGehee variables from the spherical variables with respect to the y axis.We avoid high inclinations just taking into account inclination values in [−π/4, π/4].

Figure 2 :
Figure 2: Local charts considered avoiding the z axis (left), x axis (middle) and y axis (right).The inclination considered for each local chart is always in the interval [−π/4, π/4].

Figure 3 :
Figure 3: Spherical caps not parametrized with the different local charts.Any pair of local charts recovers the whole sphere S 2 .

Figure 4 :
Figure 4: The four 1-EC orbits ((x, y) variables) that exist for the planar RTBP for values of µ = 0.1, 0.5 and C = 5 (from left to right).Two 1-EC orbits each one symmetrical with respect to the x axis (the blue ones), and the other two EC orbits are symmetrical images of one another (the red ones).

Figure 6 :
Figure 6: Left.Collision manifold Λ for the planar RTBP.Middle.Three ejection orbits, the purple, the black and the yellow ones have angular momentum m z > 0, 0, < 0 respectively at the first crossing with Σ m .The black one is an EC orbit.Right.Schematic behaviour of two orbits that go close to collision.

Figure 7 :
Figure 7: Left.Schematic plot close to collision.Middle.A set of ejection orbits integrated up to the first crossing with Σ m .The interior black curve is an EC orbit.Right.Zoom region around the origin.Set of ejection orbits and the same set close to collision.

Figure 13 :
Figure 13: µ = 0.1, C = 4 (equivalently H = −2).Left.The blue and red points (A 1 ∩ D 1 ) provide initial conditions of the two planar symmetrical and the two planar non symmetrical 1-EC orbits, respectively.The green points (A 2 ∩ D 1 ) provide initial conditions of spatial 1-EC orbits.The grey points provide initial conditions of candidates to 1-EC orbits.Right.Minimum distance at the first approach.

Figure 15 :
Figure 15: µ = 0.1.Configuration space (x, y, z).The two planar symmetrical 1-EC orbits in blue.The two non symmetrical planar ones in red.The spatial ones in green.The candidates to 1-EC orbits in grey.C = C L1 (left) and C = 6 (right).