Dynamics of the parabolic restricted three-body problem

The main purpose of the paper is the study of the motion of a massless body attracted, under the Newton’s law of gravitation, by two equal masses moving in parabolic orbits all over in the same plane, the planar parabolic restricted three-body problem. We consider the system relative to a rotating and pulsating frame where the equal masses (primaries) remain at rest. The system is gradient-like and has exactly ten hyperbolic equilibrium points lying on the boundary invariant manifolds corresponding to escape of the primaries in past and future time. The global flow of the system is described in terms of the final evolution (forwards and backwards in time) of the solutions. The invariant manifolds of the equilibrium points play a key role in the dynamics. We study the connections, restricted to the invariant boundaries, between the invariant manifolds associated to the equilibrium points. Finally we study numerically the connections in the whole phase space, paying special attention to capture and escape orbits.


Introduction
Astronomy textbooks typically present galaxies as calm, solitary and majestic island worlds of glittering stars. However, the Hubble images support the well-known idea that galaxies are dynamic and energetic. In particular, the bridges and tails seen in some multiple galaxies are just relics of close encoun- 5 ters. The consequences of the brief but violent tidal forces have been studied by Toomre and Toomre in [16] and by Namboodiri et al. in [10] considering a simple-minded fashion: each encounter is considered to involve only two galaxies describing a roughly parabolic path.
This approach of the dynamics of the close encounters for two galaxies has 10 been used, for example, by Condon et al. in [2], in the case of the galaxies UGC 12914 and 12915, or by Günthardt et al. in [7], for the system AM1003-435. The parabolic model has been also used in the study of the formation of planetary systems. Fragner and Nelson,in [6], examine the effect of parabolic stellar encounters on the evolution of a Jovian-mass giant planet forming within 15 a protoplanetary disc. Pfarzner et al., in [12], study the close encounter of two stars, one of them surrounded by a disc. More recently, in [14], Steinhausen et al. present a numerical investigation to study the effect of gravitational stardisc interactions on the disc-mass distribution, considering coplanar, prograde encounters on parabolic orbits. 20 A close approach of two galaxies (or stars surrounded by a disc) cause significant modification of the mass distribution or disc structure. Focussing just on one particle that initially stays in one galaxy (or around one star), after the close encounter, it can jump to the other galaxy or escape. One aim is to study the regions in the phase space where the particle remains or not around each 25 galaxy. To perform this study, we consider a very simple model, the so called planar parabolic restricted three-body problem, which describes the motion of a massless particle submitted to the gravitational attraction of two masses -called primaries-that move in parabolic orbits around their common center of mass, when the primaries and the particle move in the same plane. 30 As far as the authors know, very little literature has been devoted to restricted three-body problems where the gravitational attraction of the primaries is non-periodic. That is, when the energy of the primaries is non-negative, so that they move in parabolic or hyperbolic orbits. Meyer and Wang,in [9], studied restricted isosceles problems (the three bodies are at the vertices of an 35 isosceles triangle) when the energy of the primaries is non-negative and the infinitesimal mass moves in a line perpendicular to the orbital plane passing through the center of mass of the primaries. After that, Cors and Llibre, in [3] (see also the references therein), obtained a classification of the orbits in the parabolic case in terms of the asymptotic velocity of the infinitesimal mass 40 when t → ±∞ and the number of times that the infinitesimal mass intersects the plane which contains the motion of the primaries. More recently, Alvarez et al., in [1], describe some features of the planar parabolic restricted three-body problem and show the existence of special types of motion. Finally, Faintich, in [5], considers the planar hyperbolic restricted three-body problem and applies 45 the model to a hypothetical star-Sun-comet system to determine the effect of a stellar encounter on the orbit of the comet.
Being the previous works our original motivation, the main purpose of the paper is to continue and complete the results obtained in [1] and to describe the flow of the planar parabolic restricted problem (or simply parabolic problem 50 along the paper). We consider in this work the simplest case, when the two primaries have equal masses. In particular, we will focus on the role of the invariant manifolds to explain how a particle can be captured around one primary or escape.
The paper is organized as follows: in Section 2, the equations of the motion 55 of the parabolic problem are given both in an inertial system of coordinates and in a rotating and pulsating (synodic) frame. After that, in order to extend the flow of the system when the primaries are at infinity, the phase space is compactified in the time direction with the introduction of a suitable variable θ ∈ [−π/2, π/2], and the so called global system is obtained. The flow when 60 θ = ±π/2 is invariant, obtaining the upper (θ = π/2) and the lower (θ = −π/2) boundary problems. In particular the equilibrium points and the homothetic solutions associated to them are considered and the gradient-like property of the equations is remarked. Finally the possible regions of motion (the so called Hill's regions) are described as well. Section 3 is divided in two parts. First, we analyze the final evolutions for the motion of the particle. We focus on escape and capture orbits in the inertial system of coordinates and we introduce a criterium (the so called C-criterium) that allows us to classify an orbit in the synodic system of coordinates. Second, we analyze the dynamics on the boundary problems. We will see that it is crucial to know what we call map 70 of heteroclinic connections between invariant objects of the boundary system, that is, trajectories that start and end at a primary and/or an equilibrium point and/or the infinity. We recall some connections already known in [1] and we show the existence of new heteroclinic connections. Section 4 is devoted to numerical explorations. We present two strategies to find initial conditions of 75 collision orbits, both forwards and backwards in time, called connecting orbits.
Applying one or another strategy we obtain non symmetrical and symmetrical connecting orbits respectively. We will show that the equilibrium points and the invariant manifolds associated with them play a key role on the specific path of such connecting orbits and we will find orbits with close paths to triangular 80 and/or collinear configuration during their trajectory.

Description of the problem and main features
In order to have a self contained paper, we present the equations of the motion of the problem and other main features. The details can be found in [1].

85
Let us consider three bodies in an inertial (sidereal) reference system. Two of the bodies, called primaries, with masses m 1 and m 2 , move in parabolic orbits of the two body problem around their common center of mass. The third body is a particle of infinitesimal mass m 0 that moves under the gravitational attraction of the primaries without affecting them in the same plane of the motion of the 90 primaries. The problem of the description of the motion of the particle is the planar parabolic restricted three body problem (simply parabolic problem along the paper). We will consider in this paper the case of two equal masses for the primaries, that in a suitable units means that we can take m 1 = m 2 = 1/2. It is well known (see for example [4]) that, when the primaries move in a parabolic motion, the relative position vector from m 1 to m 2 is R = (σ 2 − 1, 2σ) where σ = tan(f /2), f is the true anomaly, their mutual distance is r = σ 2 + 1, 3 3 , being T the time of passage at the pericenter (see Figure 1). Then, the equation of the motion of the particle in the inertial where˙= d dt is the derivative with respect to the time t, Z 1 and Z 2 are the 95 position vectors of the primaries, and we have assumed that the constant of gravitation G = 1. After placing the common center of mass of the primaries at the origin, we perform two changes of variables. First, a standard rotating and pulsating (synodic) coordinate system z = (x, y) is introduced, so that the primaries remain fixed along the x-axis at z 1 = (− 1 2 , 0) and z 2 = ( 1 2 , 0) respectively. This change is done via the complex product (2) Second, we introduce the reparametrization of time dt ds = √ 2 r 3/2 . After some straightforward computations, the variable σ, that gives the relative position of the primaries, can be expressed in terms of the new independent variable s as σ = sinh(s), and the equations of motion for the particle in the new coordinate where ′ = d ds denote the derivative with respect to s, and Ω, the potential function, is given by Clearly, when the primaries tend to infinity along their parabolic orbits, t → ±∞ and the new time s also tends to ±∞. In order to extend the flow of the system when the primaries are at infinity, a new variable θ is introduced through the change sin(θ) = tanh(s). With the new variables (θ, z, w) the system (3) becomes the following autonomous system where Notice that the the original system (3) is defined for s ∈ (−∞, ∞), which corresponds to θ ∈ (−π/2, π/2), whereas the extended system (4) is defined also for θ = ±π/2, and it is invariant at the boundaries θ = ±π/2. Therefore the From now on, we will call system (4) the global system, and we denote as configuration space the projection of D on to the (x, y) plane. The two invariant systems will be denoted as the upper (θ = π/2) and lower (θ = −π/2) boundary problem respectively and the corresponding equations are    z ′ = w, Finally, notice that the two changes

Symmetries
The global system (4) has the following two symmetries The symmetry (7) implies that given a solution of the global system there exists another one which is symmetric with respect to y = 0 in the (x, y) plane (reversing time and θ). The symmetry (8) implies that given a solution of the 110 global system there exists another one which is symmetric with respect to the origin in the (x, y) plane.
The symmetries on the upper and lower boundary problems (5) are For each solution of a given boundary problem there exists another one of the same problem that is symmetric with respect to y = 0 in the (x, y) plane (using 115 (9)) and another one symmetric with respect to the origin (using (10)).
We emphasize that the symmetries (8) and (10) are specific for the parabolic problem with equal masses, whereas (7) and (9) apply for the parabolic problem for any value of the mass parameter.

Equilibrium points
The equilibrium points of the global system (4) are given by Therefore, all the equilibrium points of the global system are at the upper and lower boundary problems (5). Moreover, the potential function Ω satisfies whereΩ(z) is the potential function of the circular restricted three-body problem in rotating coordinates (see [15]). Thus, the equilibrium points of the parabolic problem in each boundary coincide with the classical equilibrium points of the restricted circular three-body problem: three collinear and two triangular. We denote by L + i and L − i , i = 1, ..., 5, the equilibrium points for θ = π/2, and 125 θ = −π/2 respectively. Note that due to the symmetries (9) and (10) L ± 1 and L ± 5 are opposite to L ± 3 and L ± 4 respectively. See Table 1. Linearizing the upper boundary problem, the eigenvalues associated to the collinear equilibrium points L + i , i = 1, 2, 3, are Thus, the collinear equilibrium points have an unstable manifold W u (L + i ) of dimension 1 and a stable manifold W s (L + i ) of dimension 3. Due to symmetry (6), in the lower boundary problem the dimensions are the contrary (the unstable is of dimension 3, the stable of dimension 1). In the case of the triangular equilibrium points L + i , i = 4, 5, the linearization gives the eigenvalues Thus, both invariant manifolds associated to the triangular equilibrium points, the unstable and the stable one, are of dimension 2.
Considering the equilibrium points in the global system (4), the dimension . . , 5 increases in one. See Table 1.

Homothetic solutions
Besides the equilibrium points, the simplest solutions of the global system (4) are the five homothetic solutions connecting the equilibrium points where the coordinates (x i , y i ) are given in Table 1.
Clearly these five homothetic solutions in the rotating-pulsating coordinate system are homographic solutions of the original coordinate system (1). They

Jacobi function and Hill's regions of motion
We consider the analogous function to the Jacobi constant of the circular restricted three-body problem (see [15]), that we call, by similarity, the Jacobi function: In the parabolic problem, however, C is not constant along the solutions of the global system (4), because (see Proposition 1 in [1]) Notice that C does not depend on θ explicitly, while its derivative does. Furthermore, the Jacobi function has a piecewise monotone behavior along the solutions of the global system (except at the homothetic ones, where the Jacobi function is constant). This behavior is known as gradient-like property. More precisely, along any solution of the global system, when θ ∈ [−π/2, 0] (s ≤ 0) the function C decreases, whereas for θ ∈ [0, π/2] (s ≥ 0) the function C increases. In the 145 boundary problems, the function C is monotone (for θ = −π/2, C decreases and for θ = π/2, C increases). A particular and important consequence of this gradient-like property is that there cannot exist periodic orbits.
From (12) we obtain the so called Hill's regions, that is, the allowed regions of motion in the configuration space. For a fixed value of C, we consider the zero velocity set, V 0 (C), in the configuration space defined by It can be shown that V 0 (C) is a set of closed curves (called zero velocity curves or zvc). Due to (11), the topology of the zvc is the same as in the circular restricted 150 three body problem (see [15]). In order to describe them, let C i = C(L ± i ), i = 1, . . . , 5 be the value of the Jacobi function at the equilibrium points. We observe that C 1 = C 3 , C 4 = C 5 , and in Table 1 the (approximate) values of C 2 , C 3 and C 4 are given. We plot in Figure 2 the zvc and the forbidden regions of motion (shaded regions) for different fixed values of C.

155
As the Jacobi function varies with time, the Hill's regions evolve also with time. We describe the evolution of the geometry of the zvc as the value of the Jacobi function increases. For C < C 4 = C 5 , the set V 0 (C) is empty and the motion is possible everywhere. For C 4 ≤ C < C 1,3 , the zvc are two closed curves surrounding the triangular equilibrium points (Figure 2, a) and b)), 160 which increase in size as C increases. For C 1,3 < C < C 2 , the zvc still have two components, one around the primaries and L ± 2 and an exterior one enclosing all the equilibrium points (Figure 2, c)), so there are two regions where the motion is admissible. Finally, for C > C 2 , the interior component of the zvc becomes two unconnected curves, each one surrounding one primary (see Figure 2, d)). space for fixed values of the Jacobi function. From left to right: C = 6, C = 6.9 < C 1 = C 3 , Taking into account the gradient-like property (13), the piecewise monotone behavior of the Jacobi function for s ≤ 0 and s ≥ 0, and the geometry of the zvc, we have that following a solution forwards in time, the regions of admissible motion grow when θ ∈ [−π/2, 0] and shrink when θ ∈ [0, π/2].

Dynamics of the parabolic problem 170
In order to describe the dynamics of the parabolic problem, we will focus on two aspects: the final evolutions when time tends to infinity, and the existence of heteroclinic connections between equilibrium points.
We will see that the behavior of the trajectories as time tends to infinity is rather simple, as it cannot be otherwise in a problem where periodic orbits do First, in order to have a tool to classify the orbits depending on their final evolution, we give a criterium that allows to decide when an orbit is captured by one of the primaries or when it tends to infinity. Next, we show the exis-190 tence of some heteroclinic connections and complete the diagram of heteroclinic connections given in [1].

Final evolutions
We focus on the ultimate behavior of the motion of a particle when the time tends forwards to infinity, i.e. the ω-limit (see, for example, [11]) of the We focus on escape and capture orbits. Essentially, an escape orbit will be a path along which the particle moves away from both primaries and their center 200 of mass, whereas a capture orbit will be a trajectory along which the particle approaches one of the primaries, in the sense that while the primary tends to infinity along its parabolic orbit, the particle follows the primary (spinning around it) at a bounded distance. From Painleve's theorem (see for example Chapter 4 in [13]), if the right maximal interval of existence [0, β) is finite, a collision 205 in finite time occurs. That case can be considered as capture. For this reason, from now on we will only consider orbits such that its right maximal interval of existence is [0, ∞). All the definitions and results are given considering that we can take the limit when time (t and s) tends to infinity, so that θ tends to π/2.
The definitions of capture and escape that we use in the present paper are 210 the following ones.
Definition 1. Let Z(t) be a solution of the parabolic problem given by equations (1). We say that • it is a capture orbit around the primary of mass m i , for i = 1 or 2, if It is important to point out that the definition of capture and escape is done in the inertial reference system. From the change of variables (2) we have that where the distance between the primaries r ↗ ∞ as time tends to infinity.
Clearly, to be a capture orbit it is necessary that lim s→∞ |z(s) − z i | = 0, but that condition is not sufficient. On the contrary, to be an escape orbit it is 220 enough that the lim inf s→∞ |z(s)| ≥ K, for K big enough to ensure that the particle does not approach the primaries when s ↗ ∞. The aim is to have a criterium that allows to check whether a solution satisfies these conditions in the synodic reference system.
We introduce the concept of a collision orbit in the synodic reference system. Next, our aim is to prove that in the synodic system, the solutions of the global system can be classified in three classes: collision orbits, trajectories that 230 go to infinity or trajectories tending to an equilibrium point. The first ones are candidates to be capture orbits, while the second ones are escape orbits. On one hand, θ(s k ) → π/2, so {C(γ(s k ))} k∈N is an increasing sequence. On the other hand, {Ω(z(s k )} k∈N is bounded. Thus, {|w(s k )|} k∈N must be also bounded, and there exists the limit lim k→∞ C(γ(s k )) = C ∞ < ∞.
Furthermore, {γ(s k )} k∈N is contained in a compact set, so there must exist a convergent subsequence. To simplify notation, we suppose that the sequence is convergent, so where q = (π/2, z q , w q ), belongs to the ω-limit of γ(s), and |z q | ≤ M and 235 |z q − z i | ≥ δ, ∀ i = 1, 2. Denote by Γ the solution of the global system through q.
Let B a ball of center q and radius δ/2, so that the vector field of equations As a consequence, we obtain the following results.
be a solution of the global system (4). Then, either it is a collision orbit, or lim s→∞ |z(s)| = ∞ or its ω-limit is an equilibrium point.

250
Using the above statements, we can give the following criterium, named the C-criterium because the gradient-like property of the Jacobi function C is the main key. In the second assumption, suppose that lim inf s→∞ |z(s)| ≤ K for some constant K. Again, applying the Lemma 1 we will get a contradiction. Thus, the limit lim s→∞ |z(s)| = ∞ and then it is an escape orbit.
In Theorem 4 of [1], the authors give a similar criterium (with the same name) but with some important differences. The first one is that the conclusion of the theorem is that the orbit is a collision orbit or its ω-limit is an equilibrium point, but they forgot to include the hypothesis that the particle must be in the bounded component of the Hill's region. So, they do not consider the escape orbits. Another important issue is that the authors use in the proof the fact that the Jacobi function C(γ(s)) ↗ ∞ for orbits that do not tend to an equilibrium point, although, as far as we know, this fact is not proved. Finally, the authors use the result to prove the existence of capture orbits. But we want to remark that whereas escape orbits can be detected studying their behavior in the synodic reference system, it is not sufficient to know that an orbit is a collision orbit to ensure that it is a capture orbit. In fact, in order to have a capture orbit it would be necessary to see if . On the other hand, the connections given in [1] were proved using the invariance of the planes y = y ′ = 0 and x = x ′ = 0, and the study of the flow restricted to each plane. From that, it is clear that W u (L + i ), i = 1, 2, 3, which are of dimension 1, are contained in the y = 0 axis. Then, neither L + 1 , nor L + 3 , connects with L + 2 . For the same 290 reason, there is no connection from L + 1 to the second primary or from L + 3 to the first primary. And similarly, there is no connection from L + 2 to infinity. Therefore, all possible connections starting at one of the collinear equilibrium points are already known and drawn in Figure 3. Looking at the possible connections emanating from L + 4 and L + 5 , we want to show that the following 295 heteroclinics exist: (i) from L + 4,5 to ∞ (dashed green lines), (ii) from L + 4,5 to L + 1,3 (continuous red lines), (iii) from L + 4,5 to the primaries (continuous red lines).
As we will see, we use essentially the behavior of the zero velocity sets 300 V 0 (C) and the invariant manifolds of the equilibrium points to show the above heteroclinic connections emanating from the triangular equilibrium points. A similar methodology can not be applied in the case of connections arriving at the triangular points. When exploring the behavior of the orbits of the W s (L + 4,5 ) backwards in time, the Jacobi function decreases to values lower than 4, and all 305 the configuration space is available. So, the zvc do not play any role. For that reason, we have not explored the connections from one primary, or infinity, to the triangular equilibrium points.
We start justifying connection (i): from L + 4,5 to infinity (dashed green lines in Figure 3). They can be proved immediately from the fact that the plane (y, y ′ ) is invariant, and the equilibrium points L + 4,5 are saddle. That implies, in particular, that there exist orbits of the invariant manifolds W u (L + 4,5 ) that lie on the x = 0 axis, one branch of them tending to L + 2 , the other one escaping to infinity.
Next, we want to show numerically the existence of connection (ii): from L + 4 315 to L + 3 . By symmetry, all the other connections from L + 4,5 to L + 1,3 are obtained (red continuous lines in Figure 3). We examine the behavior of the invariant manifold W u (L + 4 ) forwards in time. In order to do so, we consider an approximated parametrization of the invariant manifold, which is of dimension two, and we propagate the initial conditions given by the parametrization forwards 320 in time until a certain section close to L + 3 . To obtain such intersections with a good accuracy, the linear approximation of the parametrization is not good enough. We have considered the fourth order approximation (see [8] for details on the computation of the invariant manifold of an equilibrium point up to a desired order).

325
Given a value of the Jacobi function C = C * , we define the section and the intersection of the invariant manifold with that section is the set The invariant manifold W u (L + 4 ) is of dimension 2, so for values of C * greater but close to C 4 (on the upper boundary problem, the Jacobi function increases), the set W u C * (L + 4 ) is a simple closed curve (see Figure 4). We want to study the evolution of the sets W u C * (L + 4 ) as C * approaches C 3 ≃ 6.91359245. On one hand, the zero velocity set V 0 (C * ) has two components 330 (one in y > 0, the other one in y < 0) that approach to each other at the point Figure 2 a) and b)). On the other hand, the projection in configuration space of W u C * (L + 4 ) is a closed curve that surrounds one of the components of V 0 (C * ) (the one in the halfplane y > 0), and, when C * ↗ C 3 , this curve is trapped between the two components of V 0 (C * ). So at C * = C 3 , the 335 two components of the zero velocity curve and the set W u C3 (L + 4 ) must coincide at L + 3 , and at least one orbit on the invariant manifold must end up at the equilibrium point L + 3 (see Figure 5 and 6).  Finally, we observe that when we consider the points of the subset of W u C3 (L + 4 ) belonging to the bounded component of the Hill's region, the associated trajec-340 tories will end up at a collision with one primary except the ones that tend to L + 2 . So the heteroclinic connections (iii) from L + 4,5 to the primaries follow.  Figure 6: Upper boundary problem: projection of W u C * (L + 4 ) (blue) in configuration space and Z 0 (C * ) (black) for the values C * = 6.913 and C * = 6.9136.

Numerical results
As stated before, the simplest solutions in the parabolic problem are the homographic solutions, which remain fixed in the (x, y) plane. Among all other 345 possible orbits or trajectories in the parabolic problem, we are interested, due to astronomical reasons, in ones that go from primary to primary (either the same one or not) that is, capture orbits forwards and backwards in time. As explained in Section 3, we have a criterium to determine when a solution of the parabolic problem in the synodic system satisfies the necessary condition to be a capture 350 orbit in the sidereal system. That is, when an orbit in the synodic system is of collision forwards in time. We want to use that criterium to find sets of initial conditions corresponding to collision orbits both forwards and backwards. We call them connecting orbits. Our purpose is to present strategies that allow us to find initial conditions of such type of motion. We will also discuss strategies 355 to find connecting orbits with a priori desired path with close approaches to triangular or/and collinear configurations.
We will see that the two strategies presented are based on the knowledge of the flow on the boundary manifolds, θ = ±π/2, and on the invariant manifolds associated to their equilibrium points, and whereas the first strategy gives 360 rise to (typically) non symmetrical connecting orbits, with the second one only symmetrical connecting orbits are obtained.

Connecting orbits with passages close to collinear or triangular configurations
We consider from now on the global system given by (4) and we focus our 365 attention on connecting orbits with a passage close to a collinear or a triangular configuration. We say that a connecting orbit for the global system is of type 2} and k ∈ {1, ..., 5}, if it is a collision orbit with m i backwards in time, collision orbit with m j forwards in time, and along its trajectory it has a close approach to the homothetic solution generated by L ± k 370 -in short, we will say that the orbit has a close passage to L k , although strictly speaking there are no equilibrium points in Int(D)-. This implies that at some time the position of the particle and the primaries is close to a collinear or triangular configuration. We observe that, due to the connection diagram of each boundary problem (θ = ±π/2) discussed in the previous section, it is 375 natural to find connecting orbits that have several passages close to different L m , m ∈ {1, ..., 5}.
The procedure to find initial conditions that lead to (typically non-symmetric) connecting orbits of type m i − L k − m j is described as follows. Consider an invariant manifold of an equilibrium point of the lower boundary, W s,u (L − k ) for 380 some k ∈ {1, .., 5}. We take initial conditions (z 0 , w 0 ) on W s,u (L − k ) but with a value of θ 0 at a certain distance of −π/2, θ 0 = −π/2 + δ with δ > 0. We integrate this initial condition for the global system forwards and backwards in time for θ ∈ [−π/2 + ϵ 1 , π/2 − ϵ 2 ], with 0 < ϵ 1 , ϵ 2 , and ϵ 1 < δ. The key point is to find suitable initial conditions (θ 0 , z 0 , w 0 ) such that we can guarantee that 385 for the global flow, the associated orbit is a connecting orbit having a close path to L k (or even to more equilibrium points). In general, we take initial conditions on W s (L − k ) and a small value of δ (usually of order 10 −3 ) . Starting close to the lower boundary problem and integrating backwards, we can use the knowledge of the behavior of the invariant manifolds in that boundary problem 390 to obtain collision orbits with a desired primary. To select the appropriate initial conditions corresponding to collision orbits forwards in time, we apply the C-criterium. It is also possible to start with W u (L − i ) (that is the case of the examples shown in Figure 8). In that case, we have considered bigger values of δ (of order 0.05), so that the initial conditions are not so close to the lower 395 boundary problem. In this case is not possible to use the knowledge of the behavior of that invariant manifold on the lower boundary problem, and the C-criterium is also used backwards in time to select the initial conditions that correspond to collision orbits.
Finally, the same strategy can be applied starting at initial conditions close  . On the left, we can see the projection in the synodic configuration space, on the right, the projection in the Z coordinates. Initially, for a range of time including t 1 , the particle surrounds the 415 fixed position of m 2 in synodic coordinates (the particle spins arround m 2 on and on along its parabolic orbit in sidereal coordinates). At time t 2 the particle is close to the collinear point L 2 in synodic coordinates, which corresponds to a quasi collinear configuration where the three bodies are almost aligned. As the time increases, the particle gets captured around m 2 . (with θ close to −π/2). We want to stress that whereas along the span of time from capture around m 1 to the close passage to L 2 both orbits are quite close to each other, after the passage close to the collinear configuration, they evolve 425 in a complete different way. One orbit is captured by m 2 whereas the other one is captured by m 1 . So we remark that, some times, a small change in the initial conditions may give rise to rather different behaviors. Of course this is due to the presence of the invariant manifold of L 2 .
In Figure 9 two solutions with close passages to two equilibrium points are In this case, the orbits are almost coincident for t ≤ t 2 , leaving a neighborhood of m 2 and having a passage close to a collinear configuration. Then, they have different evolutions forwards in time. The invariant manifold of L + 2 is the responsible for the separation of both trajectories: while one particle goes to 435 collision with m 1 , the other one goes to collision with m 2 . Clearly, in between there must exists a trajectory belonging to W s (L + 2 ) connecting m 2 with L + 2 . Other examples are shown in Figure 10. The left plot shows two orbits visiting three equilibrium points: orbit "a" is a connecting orbit of type m 2 − L 2 − L 4 − L 3 − m 2 whereas orbit "b" is an escape one. It is also clear that there to −π/2, we show in Figure 11 two orbits with a close passage to a triangular

Symmetric connecting orbits
The procedure to find initial conditions that lead to symmetric connecting 450 orbits is based on symmetries (7) and (8). On one hand, using the symmetry (7) if an orbit leaves a neighborhood of one primary, and crosses the section θ = 0 such that y = x ′ = 0, the trajectory will end at the same primary. This is a connection m i − m i symmetric with respect to the x axis on the (x, y) plane.
On the other hand, using both (7) and (8), if an orbit leaves a neighborhood of 455 one primary, and crosses the section θ = 0 such that x = y ′ = 0, the trajectory will end at the other primary. This is a connection m i − m j , i ̸ = j symmetric with respect to the y axis on the (x, y) plane.
We have explored these two kind of connections in the following way. Take initial conditions at θ = 0 in the plane (x, y ′ ) (connecting orbits for the same 460 primary) and in the plane (y, x ′ ) (connecting orbits between different primaries).
Integrate such initial conditions forwards and classify them, using the C-criterium, depending whether the orbit is a collision orbit with one of the primaries or the orbit escapes to infinity.
Notice that in some cases the integration cannot be done until C = 8, when 465 the C-criterium is applied, due to the fact that the orbit has a close encounter to one primary. At these points a binary collision regularization should be performed in order to continue the integration. Instead, we classify the orbit as "possible collision with the primary".
To avoid unnecessary computations, we discard initial conditions of the 470 planes (x, y ′ ) and (y, x ′ ) at θ = 0 such that the value of the Jacobi function C is greater than C 3 or C 2 , so that the C-criterium applies directly without integration. More precisely, the curves C = C 3 and C = C 2 determine a subset of the escape and capture regions. Notice that on the plane (y, x ′ ) the curve C = C 2 does not play any role. Moreover, using the symmetries of the problem, 475 it is enough to do the computation in the positive x and y planes respectively.
The detailed regions of escape and capture on the (x, y ′ ) and (y, x ′ ) planes are shown in Figure 12, left and right, respectively.
As we saw in the previous section, the codimension 1 invariant manifolds associated to the collinear equilibrium points separate the different types of or-480 bits and these are precisely the boundaries of the escape and capture regions plotted. In order to show this assertion, we consider a region where three different regions (escape, capture around m 1 and capture around m 2 ) meet, for example in the (y, x ′ ) plane, see Figure 13 left. More precisely, we focus on regions meet, see Figure 13 right. We take a circle centered at the common point to the three different regions and a small radius, for instance 0.00005. We take initial conditions on this circle parameterized by an angle φ varying from 0 to 2π counterclockwise, and consider the corresponding orbits. We show the final evolution of these orbits varying φ in Figure 14. We start at a given point, la- beled by "a". This point belongs to the region of capture orbit around m 2 . As φ increases, the orbits are of the same type until the boundary between this region and the escape region is crossed, see orbits labeled "b" and "c". We can see in Figure 14 that there is a precise value of φ, such that the orbit tends asymptotically to L 3 . That orbit belongs to W s (L + 3 ), which separates the capture orbits 495 around m 2 and escape orbits. Increasing again φ, we obtain escape orbits until the boundary between escape and capture around m 1 is crossed, orbits "d" and "e". We can see that the point in that boundary corresponds to an orbit that belongs to W s (L + 1 ). Finally, increasing φ we have points in the capture region around m 1 until the boundary with the capture region around m 2 is reached 500 again, orbits "f" and "a". In this case, the boundary corresponds to an orbit on W s (L + 2 ). Therefore, the manifolds of L + i , i = 1, 2, 3 are the boundaries that separate the different behaviors of capture and escape. The initial conditions of the orbits labeled as "a,". . .,"f" are shown in Table 2. We also want to stress that the common point to the three regions is precisely associated to an orbit 505 tending to L + 5 . Notice that all the orbits have a close passage to L 5 . in the (y, x ′ ) plane. The regions in green and light blue correspond to points that reach a neighborhood of m 1 or m 2 , respectively, but the integration has stopped without classifying the point. The plot on the right is a magnification.

Conclusions
All possible trajectories connecting A and B (A and B being an equilibrium point, or a collision with a primary, or escape) have been established, except the ones that require a regularization of the singularities of the equations of motion.

510
This was done using the map of heteroclinic connections and the C-criterium.  Table 2: Initial conditions corresponding to the orbits a,.
manifolds associated to the equilibrium points while in the second one, the most important role is played by the symmetries of the problem.

515
In the case of interacting disk galaxies (each galaxy treated as a point mass embedded by a disk of test or massless particles) the present model shows that, after a galaxy close encounter, the regions of the phase space where the test particles remain or not around each galaxy are confined by the invariant manifolds of the equilibrium points.