HOMOCLINIC AND HETEROCLINIC TRANSFER TRAJECTORIES BETWEEN PLANAR LYAPUNOV ORBITS IN THE SUN-EARTH AND EARTH-MOON SYSTEMS

. In this paper a method for ﬁnding homoclinic and heteroclinic connections between Lyapunov orbits using invariant manifolds in a given energy surface of the planar restricted circular three body problem is developed. Moreover, the systematic application of this method to a range of Jacobi constants provides a classiﬁcation of the connections in bifurcation families. The models used correspond to the Sun-Earth+Moon and the Earth-Moon cases.


1.
Introduction. Interest in the scientific advantages of the Lagrange libration points for space missions has been increasing in the last 25 years. In 1978 ISEE-3 was launched and inserted into a nearly-periodic halo orbit around the Sun-Earth L 1 Lagrange point to pursue studies of the Earth-Sun interactions (see [1]). Afterwards, the spacecraft visited the vicinity of the L 2 libration point to explore the magnetotail of the Earth. It was then renamed as ICE (International Cometary Explorer) to have a close encounter with the comet Giacobini-Zinner. Since then, several space agencies have used libration orbits as target orbits for their probes (SOHO, MAP, ACE, Genesis. . . ) and some other future missions are also aimed at these kind of orbits (NGST, Herschel-Planck. . . ).
New space missions requirements are increasingly complex, with ever more demanding constraints and the commitment to minimize the costs. Consequently, a greater understanding of the dynamics near the collinear libration points is not less than ineluctable. The fundamental breakthrough that has given the theoretical and numerical framework for this understanding in the last decades is the use of Dynamical Systems tools. Dynamical Systems Theory, founded by Poincaré by the end of the XIX century, has used the RTBP as one of the paradigmatic models for its application, looking at this problem from a global point of view. The goal is to get a picture of the evolution of all the states of the system, looking at all the orbits together instead of individually. This approach has successfully been used in some missions such as Genesis (see [2]).
The introduction of invariant manifolds as a means to describe the phase space yields not only a more efficient determination of the desired transfers but also an adaptable procedure for mission analysis. This means that the baseline trajectory Space missions are becoming increasingly required, so knowledge about the delicate dynamics of the three body problem is essential. The spatial restricted three body problem is more difficult to tackle, due to its higher number of degrees of freedom. However, there is some evidence that starting by planar connections is a good and necessary step to be taken before moving to the 3D case (see [13]). Besides, many missions between libration points will have a small out-of-plane component; this is the case of the Earth-Moon problem. Finally, for higher inclinations the patterns of the xy projections of 3-dimensional connections are similar to the planar homoclinic and heteroclinic connections. All these facts justify that a good knowledge of the planar case is essential.
2. Methodology. The study of the motion of an infinitesimal particle affected by the gravitational attraction of two massive bodies (primaries), which describe circular orbits around their common center of masses, is called the Restricted Three Body Problem. If we only study the movement of the particle in the plane determined by the two-body movement of the primaries, we say that the problem is planar (PRTBP).
Taking adequate units of mass, length and time we can simplify the equations of motion. Furthermore, the following synodical system of reference is taken: the origin set at the center of masses of the two primaries, the X axis given by the line that goes from the smallest to the biggest primary (with this orientation) and the Y axis taken so that XY is a planar positively oriented coordinate system. In our synodical system the primary of mass µ (small) is always located at (µ − 1,0) and the primary of mass 1 − µ at (µ,0). The equations of motion are: The PRTBP in synodical coordinates has five equilibrium points. Three of them are located on the X axis, and are called Eulerian or collinear points (L 1,2,3 ). L 4 and L 5 (Lagrangian points) are found as the third vertex of the two equilateral triangles that can be formed using the primaries as vertices (figure 1).
In this paper we are interested in the dynamics concerning the equilibrium points L 1 and L 2 . The linear behavior of these points is of the type center × center × saddle. This hyperbolicity is inherited by the libration orbits, and all of them are unstable in a big neighborhood of these points. The instability can be skipped, and some bounded orbits computed, by reducing the equations to the center manifold (see [14]).
One can compute the invariant manifold structures associated to the collinear libration points and periodic orbits around them using dynamical systems theory. In particular, the stable and unstable invariant manifold tubes are asymptotic to the Lyapunov orbits. These tubes provide the key to study the natural dynamics of the libration regions. The purpose of this paper is to present the main homoclinic and heteroclinic connections between the L 1 and L 2 planar Lyapunov orbits from the point of view of practical astrodynamical applications. Orbits lingering unnecessarily about the libration points or about the small primary will be skipped.
2.1. Lindstedt-Poincaré expansions. According to Poincaré theorem, at a given energy level, there is a unique planar Lyapunov orbit homeomorphic to S 1 around each libration point (L 1 and L 2 ). This is the unique periodic motion around the L i for the planar case and the one we will be working with in order to find the connections.
The planar Lyapunov orbits and their hyperbolic manifolds can be computed using Lindstedt-Poincaré procedures. In this way one obtains their expansions in convenient RTBP coordinates.
To this end, we set the origin of coordinates at the libration point (L 1 or L 2 ) and scale the variables in such a way that the distance from the equilibrium point to the small primary is equal to one. The expansion of the equations of motion (1) in these variables (x,ȳ) takes the form, where ρ 2 =x 2 +ȳ 2 , P n is the Legendre polynomial of degree n, and c n are constants which depend only on µ and the selected equilibrium point (see [7]). We note that in (2) the linear terms appear in the left hand side part of the equations and the nonlinear ones in the right hand side. The solution of the linear part of the equations (2) is:x where κ 1 , κ 2 , ω p and λ 0 are constants for a given model and libration point. The α's are free amplitudes. α 1 and α 2 are the ones associated with the hyperbolic manifolds. If α 1 = α 2 = 0, we have the linear part of the Lyapunov orbit with amplitude α 3 . When α 1 = 0 and α 2 = 0 we have orbits tending to the Lyapunov orbit of amplitude α 3 when time tends to infinity (stable manifold). On the contrary when α 2 = 0 and α 1 = 0, orbits leave the vicinity of the Lyapunov exponentially fast in forward time (unstable manifold).
When we also consider the non-linear terms of (2), solutions are obtained by means of formal series in powers of the amplitudes of the form: where θ 1 = ωt + φ, θ 2 = λt and, Summation is extended over all i, j, k and p ∈ N. However, due to symmetries, many of the coefficients x p ijk ,x p ijk , y p ijk ,ȳ p ijk , ω ijk , λ ijk are zero. Moreover the series are truncated at a certain (high) order (see [7] for details). Nevertheless, we note that the meaning of the amplitudes in the nonlinear expansions (4) is the same one as in the linear solutions (3).

2.2.
Fixed energy surfaces. Introducing momenta p x =ẋ − y and p y =ẏ + x, the RTBP equations of motion cast into a Hamiltonian system. The Hamiltonian function is, In this paper, however, we will not be using H but the Jacobi constant, C, defined as follows, C(x, y,ẋ,ẏ) = −(ẋ 2 +ẏ 2 ) + 2 Ω(x, y) (5) with Ω as in (1). It is easily proved that C = −2 H.
The number of degrees of freedom in the planar restricted three body problem is n = 2. Thus, the order of the system is 2n = 4, which is reducible by the Jacobi constant to 3. That is to say that as C does not vary along the solutions, we have to study only 3 of the four coordinates of the phase-space for each orbit, obtaining the fourth one from equation (5) if necessary. Besides, the three dimensional restricted 3-body problem is of order 6, and there is still only one integral, so the surface of section becomes four dimensional.
The level surfaces of the Jacobi constant in the planar problem are usually called energy surfaces, three dimensional manifolds implicitly defined by equation (5) on which we can look for solutions of the PRTBP: The range of values for the Jacobi constant, [C min ,C max ], will depend on the type of connections we are looking for. When C=C Li the zero velocity curves collapse at L i , so for all bigger Jacobi constants there is no room for Lyapunov orbits around L i . Thus, we define C i max = C Li . High values of C correspond to tiny y-amplitudes of the Lyapunov orbits, as the zero velocity curves progressively become closer to the libration point.
The values C min we use are found when the Lindstedt-Poincaré series (4) no longer give accurate results when truncated at order 15. If we wanted to work with smaller C, we would have to use numerical continuation techniques. However, orbits with C near our C min are already too big to have a practical interest, with y-amplitudes around 1.1 × 10 6 km in the Sun-Earth case, and 50, 000 km in the Earth-Moon case, which are huge orbits in either problem.
To sum up, the range of Jacobi constants we work with is small, but it allows the y-amplitudes of the Lyapunov orbits to experiment a considerable variation (from almost zero up to a million km). See figure 2.2, where Lyapunov orbits for all usable C are depicted. As the Jacobi constant increases, the amplitude decreases and the other way round.

KS-Regularization.
The distance from the infinitesimal mass to both the primary bodies appears as a denominator in the RTBP equations. Approaching the biggest one is not a problem to be taken into account if we work in L 1 and L 2 regions. However, the distance to the small m i can become very short at some  points, with the consequent computational problems due to this singularity. We have used a well known regularization method introduced by Levi-Civita (2-D) and Kustanheimo in order to overcome this drawback (see [17], [18]). By using this regularization, and applying it to the integration of the RTBP vector field, we enhance the global adaptation of the model to the problem. This is to say, adding to the set of solution orbits which travel close to the Earth. These orbits, indeed, might as well be the most important in terms of practical applications for spatial missions; for instance: injection, phasing loops or reentry.
2.4. Homoclinic and heteroclinic phenomena. As stated in section 2.1, for a given energy level there is a unique planar Lyapunov orbit homeomorphic to S 1 , around each libration point (L 1 and L 2 ). As the phase space near these points has a saddle component, there are orbits asymptotically approaching the Lyapunov in forward time (stable manifolds) and orbits leaving it as well (unstable manifold). These manifolds are two dimensional in M(µ, C) (for details see [14]).
By matching one orbit from the stable manifold with another one from the unstable one, a zero cost trajectory is obtained which asymptotically approaches a Lyapunov orbit both in forward and in backward time. If the manifolds belong to the same orbit, the connections are homoclinic. On the contrary, if they come from different orbits, they are called heteroclinic connections. In our case, as there is a unique periodic orbit around each libration point for any fixed energy level, the classification in homoclinic and heteroclinic can also be regarded as trajectories involving only one L i (homoclinic) or both of them (heteroclinic).
These connections can be found by means of a Poincaré section representation. For a given M(µ, C), a Poincaré section S, can be placed at x = µ − 1 and a positive integer, k, can be chosen. We can integrate initial conditions which lay in the invariant stable and unstable manifolds of the corresponding Lyapunov orbit until the trajectory has crossed the section k times. After these crossings, we have the representation of the k-cut of the manifold with the section.
Lyapunov planar orbits are S 1 like objects. The manifolds which arise from them are like tubes in the phase space. So, they result in curves when intersected with a transversal section, S. We will write W (u/s),j i standing for the j-th intersection of the W u/s (unstable or stable invariant manifold) of the Lyapunov orbit around L i with S. Depending on the initial phase and the Jacobi constant, some orbits escape to the exterior region (see figure 3) or collide with the Earth or Moon after a maximum number of cuts. So they never reach the section again, or they do it breaking the structure of S 1 . Not all W (u/s),j i , especially as j increases, will be like S 1 (see [3]). We must state that in this paper we refer to collision when an orbit approaches the point-mass Earth or Moon within a distance slightly bigger than their respective real radii.  Once we have W (u/s),j i on S, it is convenient to look at it as a curve in the (y,ẏ) plane. In S, x is fixed andẋ can be computed using equation (5). In fact, |ẋ| is determined by the Jacobi constant, but not its sign. We need to keep in mind which is the direction of the manifolds we intersect in order to avoid problems when defining the mentioned sign across S. Thus, a point (y 0 ,ẏ 0 ) is enough to determine an initial condition once we know how to choose the sign for theẋ component. Moreover, as the coordinates we are working with give the system an autonomous character, each point in the phase space (x, y,ẋ,ẏ) determines one and only one orbit or solution to equations (1).
There exist different kinds of connections: If a point (y 0 ,ẏ 0 ) belongs to one of these intersections (see figure 4), we can complete it by finding (x 0 ,ẋ 0 ) in the Poincaré section using the aforementioned dependence on C. By unicity of the solutions of the system, this (x 0 , y 0 ,ẋ 0 ,ẏ 0 ) provides us with an orbit that leaves the vicinity of one Lyapunov orbit to approach another one. We obtain such orbit by integrating both forward and backward the intersecting point on the invariant unstable and stable manifolds respectively.
A natural way of classifying the connections consists of counting how many times they go around the small primary, the Earth or the Moon in our case. The number of loops, as a function of the number of cuts of each manifold with the Poincaré section, is (j 1 + j 2 − 1)/2 for the homoclinic orbits and (j 1 + j 2 − 2)/2 for the heteroclinic ones. Consequently, this provides a parity criterion for the total number of times the Poincaré section is crossed by every connection: So it has to be even for homoclinic trajectories and odd for heteroclinic ones.
According to the number of loops, we will note: Notation. Ho n i (n-homoclinic orbits). Homoclinic trajectory of a Lyapunov orbit around L i , with a total number of loops around the Earth equal to n (all of them travel around the Earth in the counterclockwise direction).
Notation. He n i 1 ,i 2 . Heteroclinic trajectory from a Lyapunov orbit around L i 1 to a Lyapunov orbit around L i2 winding around the small primary n times. This includes all heteroclinic connections obtained as 2.5. Some details on the numerical methodology. The numerical process for finding homoclinic and heteroclinic connections can be split into two parts: . All the implementation has been coded in FORTRAN and run in a PC with Debian-Linux operating system. 2.5.1. Integration of the manifolds. Lyapunov orbits are described by equation (4) with α 1 = α 2 = 0. The amplitude, α 3 , depends on C (for each energy surface there is one and only one α 3 , as there is a unique Lyapunov periodic orbit). In order to take initial conditions on the unstable (respectively stable) manifolds we set α 1 = (respectively α 2 = ). is a small parameter whose sign indicates towards which branch of the manifold we integrate. The phase θ 1 can be interpreted as the parametrization of the Lyapunov. On the other hand, θ 2 describes the rotation of the trajectory on the asymptotic manifold with respect to the initial position represented by θ 2 (0) = 0.
For example, if we want to find W u,2 1 for C = C * we have to set α 1 = + , α 2 = 0 and α * 3 = α 3 (C * ). Finally, we discretize the set of initial conditions by taking n θ 1 phases.
The bunch of orbits launched by these kind of initial conditions form manifold tubes like the ones represented in figure 7. A Runge Kutta Fehlberg 7-8 integrator is used. As the Poincaré section is {x = µ − 1}, the number of cuts with the section can be easily controlled by studying the sign of x(t) + 1 − µ, and refining the N-cut.

2.5.2.
Computation of the intersections between the cuts of asymptotic manifolds in the Poincaré section. Given two different manifold cuts in the Poincaré section, W u,j 1 i1 and W s,j 2 i2 , we think of them as curves (see figure 4). However, we know by the discretization taken in the initial conditions that they are polygons with n+1 sides (segments). Therefore, we can easily check for intersections between segments belonging to one of the manifold cuts and segments which form the sides of the other manifold cut.
When an intersection between the aforementioned segments has been detected, it is refined in the phase parameter (θ 1 ) of the Lyapunovs up to the required order of precision. At the end of the refining process, we store the value of the resulting phases as the representation of the connecting trajectory. The connection can then be explicitly found by integrating in forward (respectively backward) time the initial condition represented by the unstable (resp. stable) intersecting phase until the desired number of cuts with the section is reached (see figures 6, 9, 10 and 12).
3. Families of connections. One can start looking for homoclinic and heteroclinic trajectories using the elements described in the previous section. Due to the continuous dependence of the solutions with respect to the initial conditions and the asymptotic character of the manifolds, once a connection has been found there are infinitely many others in a vicinity of it (see [3]). Both by slowly varying the initial phases on the Lyapunovs or maintaining the phases and varying the Jacobi constant, similar connections can be found. In addition, a tiny modification in the initial conditions can also lead to an almost identical orbit in the phase space, which differs from the initial one in the time they spend winding around the original Lyapunov or the final one during its mid-course. Obviously, the time span for an asymptotic orbit laying in an hyperbolic manifold to leave and reach the limiting sets is infinite. So, for practical applications, by the time spent winding around such an orbit we mean the time interval in which the asymptotic connection and the Lyapunov orbit can still be distinguished as different objects at a reasonable scale.
To sum up, infinitely many connections can be found from any initial intersection on the Poincaré section. However, the classification that has been done in the following sections aims at astrodynamical applications and it deals with simple paths, in the sense that they only wind around Lyapunov orbits in the departure and arrival parts but not during the mid-course.
Furthermore, connections with a small number of loops around the Earth are usually preferred for practical applications. For instance when the intention is to move from one libration point to the other in a fast way. Nevertheless, we include homoclinic connections up to 8 loops for the Sun-Earth case, and heteroclinic connections up to 5 loops. These kinds of trajectories can be useful for a scientific observation mission, aimed at spending some time around the Earth for example.
Finally, it is also important to note that homoclinic connections cannot be found without going around the small primary at least once. This is a consequence that the manifolds never cross the x-axis before the x = µ − 1 plane (see figure 7). On the contrary, the simplest heteroclinic connections are found by intersecting the first cut of the manifolds from both sides, giving a number of loops n = 0.
3.1. Homoclinic connecting orbits. We start by C i min and find the points on S which represent a homoclinic connection on this energy surface. We can do the same for slowly increasing values of C until we reach C i max (i=1,2). Results are shown in figure 5. The y-coordinate in the Poincaré section {x = µ − 1} of the connecting trajectories is represented for every Jacobi constant. In these representations, the number of cuts of the manifolds W s , W u with S is chosen so that |j s − j u | = 1 and y > 0.
In the aforementioned figures each family of orbits is depicted using the same line style, which correspond to orbits with the same number of loops around the small primary. Each one of the families has different branches that come close to each other as C increases and finally meet at a bifurcation value, C bif . If we pick a C such that = C bif − C is big, the corresponding orbits in each branch are quite different. However, as decreases, the families approach, tending to a common limiting orbit associated with C bif . The evolution of the two branches of homoclinic connecting trajectories around L 1 in the Sun-Earth problem, Ho 1 1 , which meet at a bifurcation value for the Jacobi constant about to 3.00088389, is shown in figure 6.    On the right picture the (x, y) projection of the stable (bottom) and unstable (top) manifolds for C=3.0009 (the biggest one used for L 1 ) is shown. Again, they do not cut before crossing the Poincaré section at least once.

Symmetries in the homoclinic families. In the restricted planar three body problem, if a curve (x(t), y(t),ẋ(t),ẏ(t)) is a solution of the equations then
is also a solution. Some asymptotic orbits are a closed set with respect to this symmetry property, while some others are not. This motivates the following definitions.
Definition 2. γ 1 , γ 2 are complementary orbits (or families of orbits) if they satisfy, We note that Lyapunov periodic orbits are symmetric. Also, all 1-homoclinic orbits are symmetric as defined in point 1 above. For n-homoclinic orbits, n > 1, we can see in figure 5 that there are 4 branches in each family: two of them are symmetric (definition 1), while the other two are complementary to each other (definition 2). See figure 8 for the classification of these branches according to their symmetry properties. Around L 2 the classification of homoclinic orbits into symmetric and complementary orbits is qualitatively the same.
In figures 9 1.a and 1.b a representation of a symmetric 1-homoclinic trajectory and a pair of complementary 2-homoclinic ones for the Sun-Earth case can be found. Respectively, figures 10 1.a and 1.b show a symmetric 2-homoclinic orbit around L 1 and a pair of 3-homoclinic complementary ones around L 2 for the Earth-Moon case.
On the other hand, the planar three body problem with 0 < µ < 1 2 is not symmetric in L 1 -L 2 . However, as µ → 0 we have an increasing "almost-symmetry" with respect to the x = µ − 1 axis. The limiting case, Hill's problem, has exact geometrical symmetry with respect to this vertical axis crossing the primary. When γ 1 is homoclinic around L 1 , then its almost-vertical-symmetric partner γ 2 has to be an homoclinic trajectory around L 2 . To find which families are almost-symmetric to each other with respect to the surface of section x = µ − 1 we just have to compare figures 5.a and 5.b. If we choose a number of loops around the Earth, n, and pick one of the four branches which represent n-homoclinic trajectories around L i 1 (two if n = 1) the corresponding branch of trajectories around L i 2 contains n-homoclinic orbits which are vertically almost-symmetric to the first ones. For example, see trajectories in figure 9 1.c.
If µ is not very small vertical symmetry is lost. For instance in the Earth-Moon case. See figure 10 1.c. 3.2. Heteroclinic connections. If we intersect the hyperbolic manifolds of different L i periodic orbits, what we have is a set of heteroclinic connections instead of homoclinic orbits. Nevertheless, for all the rest, the procedure is exactly the same as before.
Obviously, the range of values of the Jacobi constant for which we can search for heteroclinic connections between Lyapunov orbits has to be the intersection between the intervals in which we have restricted our expansions for both libration points. That is [C 1 min , C 2 max ]. If we proceed exactly as we did for homoclinic orbits, storing the connections for each value of C, what we get is represented in figure 11. In more detail, the points of any vertical line (corresponding to a particular value of C) which intersect one of the curves of the figure, represent heteroclinic connections between L 1 and L 2 going through the Poincaré section with the y value represented in the y-axis of picture (i.e. y-axis corresponds to the y coordinate of the connections when x = −1 + µ). The figure is symmetric with respect to y = 0, but connections with y > 0 go from L 1 to L 2 , while connections with y < 0 go from L 2 to L 1 and   are symmetric to the previous ones due to RTBP symmetry which will be further explained for heteroclinic trajectories in section 3.2.1. The bifurcation phenomena is similar, with pairs of families of connections tending to a single one at a bifurcation value of the Jacobi constant. In figure 12 we can observe the evolution of a family of heteroclinic connections from L 1 to L 2 for the Sun-Earth problem, which has two branches that end at the bifurcation trajectory with C = 3.000863625.
We say that a pair of heteroclinic orbits, one from L 1 to L 2 and the other from L 2 to L 1 is a heteroclinic channel. Figures 9 and 10 contain some representations of heteroclinic channels.
3.2.1. Symmetries in the heteroclinic families. The intrinsic symmetry property of the PRTBP is responsible for the existence of what we have called symmetric heteroclinic channels. If (x(t), y(t),ẋ(t),ẏ(t), t) is a heteroclinic connection from L 1 to L 2 , then (x(t),ỹ(t),ẋ(t),ẏ(t)) = (x(−t), −y(−t), −ẋ(−t),ẏ(−t)), is a connection from L 2 to L 1 (see figure 11). That is, The vertical almost symmetry that exists for small values of µ does not give rise to "return" heteroclinic channels. That is to say, if γ 1 is a heteroclinic connection from L i 1 to L i 2 , and γ 2 is its vertical-almost-symmetric partner, then γ 2 also goes from L i 1 to L i 2 (figure 9 2.c). In figure 9 some examples of symmetric heteroclinic channels, as well as vertical almost-symmetric pairs of connections, are shown. For the Earth-Moon case, the almost symmetry is lost as we have seen for homoclinic connections, especially for those values of the Jacobi constant which produce Lyapunov orbits of considerably different size in each equilibrium point. If Lyapunov orbits are similar (small values of the Jacobi constant), we can find some almost vertical heteroclinic connections, as shown in figure 10 2.c.

3.3.
Homoclinic and heteroclinic connections as separatrices. For the values of the Jacobi constant considered in our studies, the zero velocity curves of the PRTBP define three regions: the region surrounding the big primary, the one surrounding the small primary and the region between L 1 and L 2 (see figure 3). Transit orbits are those which go from the region of one primary to the other. Non-transit orbits, on the contrary, bounce back to their region of origin when approaching the bottleneck of the zero velocity curve.   On the other hand, the energy surfaces of the PRTBP are 3-dimensional. Thus, invariant manifolds, which are 2-dimensional, act as separatrices of the two different types of motion: transit and non-transit orbits (see [5], [6], [9]).
Let us define κ(x, y,ẋ,ẏ) as the maximum number of cuts with the Poincaré section that the trajectory starting at (x, y,ẋ,ẏ) for t = 0 can stand, without escaping from the L 1 -L 2 region or colliding with the small primary. So far, we have searched for connections by finding the intersections between asymptotic invariant manifolds W u,j1 i 1 ∩ W s,j2 i 2 on S. These kind of banana-like plots (see figure 4) clearly divide the points (y,ẏ) in two groups: where • W denotes the interior of the curve W in S. In the first case, the points are inside the tube connecting different regions. Therefore, they escape from the L 1 -L 2 after a certain number of cuts with S. On the contrary, points which are not in the intersection of the interior of the tubes, are confined to the libration point regions and they can result in a greater number of cuts with the Poincaré section. An orbit cannot change from one regime to another, neither backwards, nor forwards in time. What can happen, however, is that it gets captured by the unstable manifolds and winds forever towards the corresponding Lyapunov orbit. If a slight change in the energy level occurs, the orbits no longer have the same regime. So, an orbit like the one we considered before, which ended up by belonging to one of the manifolds, can now belong to the transit trajectories. Consequently, κ(x, y,ẋ,ẏ) can be a good indicator to help us classify all possible initial conditions. Moreover, the homoclinic and heteroclinic connections lay in the border of In this way, they can be regarded as separatrices. For instance, let us consider the invariant unstable manifold of the Lyapunov orbit around L 1 for the Sun-Earth case which goes to the L 1 -L 2 region, for all Jacobi constants in the range with which we have been working. Using (4) with α 1 = 10 −3 , α 2 = 0 and moving φ ∈ [0, 2π] we compute κ(x, y,ẋ,ẏ). The results are represented in figure 13. In particular, the x−axis represents the Jacobi constants, and the y−axis the phase around the Lyapunov. Each point in the picture (C, φ) belongs to a coloured region whose colour represents the maximum number of cuts with the Poincaré section were possible for the particular initial condition defined by (0, 10 −3 , α 3 (C), φ). Figure 13. For the range of Jacobi constants represented in the xaxis, we consider the unstable manifold of the corresponding Lyapunov orbit around L 1 which goes to the L 1 -L 2 region. Manifolds are parametrized by the angle φ (y-axis). In the figure we represent the maximum possible number of cuts with the Poincaré section. Homoclinic and heteroclinic connections are also depicted, coinciding with the borders of the colored regions.
Once we have the coloured regions as a result of the maximum number of cuts, we can represent the connections we had previously found and depicted in figures 11 and 5, but using the y-axis to represent the phase around the Lyapunov which gives  Figure 14. Classification of the initial conditions with respect to their final behavior. Collision with the Earth (distance less than Earth radius) or escape. White zones correspond to orbits which continue to cross the section even after more than 30 cuts.
rise to the connection instead of the y−coordinate of this connection on the Poincaré section. By adding these two pieces of information (number of cuts + previously found connections), we obtain figure 13, and we can clearly see how the heteroclinic and homoclinic connections act as separatrices of the regions with constant κ. Note again that the families of connections depicted in the figure are the same as in figure  5 a.1 and 11 from L 1 to L 2 . However, they are now characterized by their departing phase in the corresponding level of energy, instead of the vertical coordinate on the Poincaré section.
The borders of the regions with an odd κ correspond to heteroclinic connections. Particularly, to the heteroclinic families which perform (κ − 1)/2 loops around the small primary. On the other hand, the homoclinic families wrap up the regions with even κ. A homoclinic connection in the border of a region with maximum number of cuts κ performs κ/2 loops around the small primary.
The points in the figure which are marked as having 0 cuts with the section correspond to trajectories which collide with the Earth before having gone through S even once. Furthermore, trajectories which are marked with a positive number of cuts can either follow the collision manifolds and get captured, or escape to another region, after having completed all the intersections with the Poincaré section. See figure 14. 4. Conclusions. In this paper, a method for numerically computing homoclinic and heteroclinic trajectories between Lyapunov orbits for any given Jacobi constant is developed. Furthermore, it is applied to a useful range of Jacobi constants, to find families of connections for the Sun-Earth and the Earth-Moon cases. The representation of these families in bifurcation diagrams provides an easy way to choose an adequate connection, in terms of requirements such as energy level or number of loops. Important symmetries are also stated and used in the computations, specially for heteroclinic connections.