Tails and bridges in the parabolic restricted three-body problem

After a close encounter of two galaxies, bridges and tails can be seen between or around them. A bridge would be a spiral arm between a galaxy and its companion, whereas a tail would correspond to a long and curving set of debris escaping from the galaxy. The goal of this paper is to present a mechanism, applying techniques of dynamical systems theory, that explains the formation of tails and bridges between galaxies in a simple model, the so-called parabolic restricted three-body problem, i.e. we study the motion of a particle under the gravitational influence of two primaries describing parabolic orbits. The equilibrium points and the final evolutions in this problem are recalled,and we show that the invariant manifolds of the collinear equilibrium points and the ones of the collision manifold explain the formation of bridges and tails. Massive numerical simulations are carried out and their application to recover previous results are also analysed.


Introduction
Gaia data release 1 has reported very recently the discovery of tails around the Large and Small Magellanic Clouds (a pair of massive dwarf galaxies) as well as an almost continuous stellar bridge between them (see [5]). Actually, in the seventies, the observation of tails and bridges in multiple galaxies was already recorded. We mention the interacting pairs M51 and NGC 5195 or the pair of interconnected galaxies Arp 295 as two particular examples (see [18] and references therein). These papers argue that tails and bridges are just tidal relics of close encounters between two galaxies. In order to study the effects of the brief but violent tidal forces due to a close encounter between the galaxies, several authors have considered a very simple model: each encounter involves only two galaxies assumed to describe parabolic orbits, and each galaxy is idealized as just a disk of non-interacting test particles which initially orbit a central mass point. This model corresponds to the parabolic restricted three-body problem (the parabolic problem along the paper), assuming that the two point primaries are the galaxies describing parabolic orbits around their common centre of mass.
There are several studies of the observable bridges and tails in galaxies. For instance, in [6] show that galaxies UGC 12914 and UGC 12915 have a continuum bridge, that is thought to be due to the collision of the galaxies 2×10 7 years ago, considering that the orbits are nearly parabolic. In [10], the authors consider the system AM1003-435, that is composed by two interacting galaxies. They studied the dynamical evolution of the encounter between the galaxies to conclude that they were moving in parabolic orbits. The N -body simulation of the orbits of stars in the galaxies shows bridges and tails. Also using the parabolic model, [12] studied the existence of bridges and tails in interacting galaxies depending on the circular velocity of the stars within the galaxies.
The parabolic model has also been used in the study of close encounters between disc-surrounded stars and the formation of planets. [13] studied the change of mass between stars when one or both of them are surrounded by a disc of low-mass particles. They concluded that, in the coplanar case, there were more change of particles between stars when the encounter was prograde. [9] studied the effect of parabolic encounters in the formation of Jovian-mass planets. They concluded that planets that have been formed after encounters are more massive and also have greater semi-major axes. [16] studied the influence on the initial density of the particles in the change of mass between star-disk encounters, concluding that the shape of the mass distribution has a high effect on the final outcome. Finally, [8], considered a Sun-star-comet system to determine the effect of the stellar encounter on the trajectory of the comet, but considering a hyperbolic model instead of a parabolic one.
The goal of this paper is, applying techniques of dynamical systems theory, to describe a mechanism that explains the formation of bridges and tails in the very simple model of the parabolic problem. Without trying to make a definition, a bridge would be an arm between the two galaxies, whereas a tail would correspond to a long and curving set of particles escaping from a galaxy.
More precisely, along the paper, we show that the invariant manifolds of the equilibrium points of the parabolic problem, and those of the equilibrium points inside the collision manifolds are the clue to find out such mechanism.
We do massive numerical simulations, considering both equal and unequal primaries, and we show the unambiguous appearance of bridges and tails, after the close encounter of the primaries. In particular, inspired in [18], we repeat some of their computations and conclude that our mechanism applies to their simulations (although in that paper, no mention of dynamical systems tools is done at all).
The application of these tools is well known in Celestial Mechanics. However, as far as we know, the application of invariant manifolds to study the close approach of two galaxies is a novelty. It is worth mentioning that there exist several papers where invariant manifolds have been applied in galactic dynamics.
In a germinal paper [15], the authors introduce a new theory for the formation of ring structures in barred galaxies. Concretely, they propose that rings are formed by material from the invariant manifolds associated to a type of periodic orbits. An study of building blocks and the morphology of rings and spirals in barred galaxies can be found in [2] and [3]. See also [19] and the references therein.
In [4] we studied the main features of the parabolic problem only for the case of equal primaries. In the present paper, we consider the parabolic problem for any value of the mass parameter µ, where 1 − µ and µ, for µ ∈ (0, 1/2], are the masses of the two primaries in normalized units. First, we show the main features of the problem and second we show the mechanism that explains the formation of tails and bridges.
More concretely, in Section 2 we describe the parabolic problem and the main relevant properties. The parabolic problem is gradient-like due to the existence of a piecewise monotone function, called Jacobi function. This property allows to classify all possible final evolutions on the dynamics of the parabolic problem.
On one hand, the flow of the system is extended when the primaries are at infinity, so the phase space is compactified in the time variable and we obtain what we call the global system. The equilibrium points of the global system and their invariant manifolds will be some of the main actors in the description of the dynamics of the problem. In particular, we show that the invariant manifolds of co-dimension 1 play a key role because they separate the different types of orbits, that is, they act as frontiers and divide the phase space in regions where the final evolution is either capture or escape. On the other hand, since we are interested in solutions that have close paths to the primaries (or even collide with them), the regularization of the equations in synodical coordinates is also performed (as far as the authors know, no paper dealing with the parabolic problem has ever considered regularized equations). There are two collision manifolds that correspond to a collision between the particle and each primary.
The associated stable/unstable invariant manifolds will be the remaining cast of actors in the paper. Section 3 is devoted to show that the stable invariant manifolds associated with the collinear equilibrium points and the unstable manifold associated to the equilibrium points in the collision manifolds are responsible for the existence of bridges and tails. The results are exemplified by some numerical explorations.
In the discussion section, some results from [18] are recovered and compared with the ones obtained by the authors. Finally we draw some conclusions.

Description of the model and main features
In order to describe the dynamical mechanisms that can help to understand the formation of bridges and tails when two galaxies have a close encounter, we consider a model already used by Toomre and Toomre. On one hand, each galaxy is modelled as a disc of non-interacting particles orbiting around a central mass point, whereas the two mass points move in parabolic orbits with respect to their centre of mass. Therefore, our simulations are based on a parabolic restricted three-body problem with a set of massless particles. The model is rather simple and does not pretend to simulate the complicated inner dynamics of the galaxies. But although its simplicity, this academic model is good enough to give a dynamical explanation of the formation of bridges and tails, that can be used as a basis to understand more complicated models.
In this Section we present briefly the equations of the motion of the model, called parabolic problem, and other main features. The details can be found in [1] and [4].

Equations of motion
Let m 1 and m 2 be two bodies, called primaries, moving in parabolic orbits around their common centre of mass. Consider a third body with infinitesimal mass moving under the gravitational attraction of the primaries in the same plane of the motion without affecting them. The planar parabolic restricted three body problem (simply parabolic problem along the paper) describes the motion of the infinitesimal mass.
We can consider, without loss of generality, suitable units of mass, length and time, such that the constant of gravitation is equal to one and, the masses of the primaries are m 1 = 1 − µ and m 2 = µ, µ ∈ (0, 0.5], called the mass parameter. Without any other interaction, when the primaries move in a parabolic motion, the relative position vector from m 2 to m 1 can be written as where the angle f is known as the true anomaly, and it is measured in the direction in which the relative position is described, starting from pericentre (see, for example, [7]). Notice that with the normalized units, the minimum distance between the primaries (when f = 0) is equal to 1. In this case, the Then, the equation of the motion of an infinitesimal mass moving under the gravitational attraction of the two point masses in an inertial system of coordinates Z = (X, Y ), with origin located at the centre of mass of the primaries, is given byZ where . = d/dt, Z 1 = Z 1 (t) = µR and Z 2 = Z 2 (t) = (µ − 1)R, are the position vectors of the primaries, see Figure 1.
First, we introduce a rotating and pulsating coordinate system z = (x, y) where the primaries remain fixed along the new x-axis at z 1 = (µ, 0) and z 2 = (µ − 1, 0). The change is given by the complex product Second, a change of time via dt ds = √ 2 R 3/2 , where R = |R| is introduced. The variable σ can be expressed in terms of the new time s as σ = sinh(s). Notice that the primaries tend to infinity along their parabolic orbits, that is, R → ∞ when t → ±∞, or equivalently, when s tends to ±∞, and that the closest passage takes place at s = 0. The system obtained with such changes is nonautonomous and non-periodic, that is, the vector field depends on non-periodic functions, concretely, on tanh(s) and sech(s).
In order to apply the standard tools of dynamical systems, like the study of the existence of equilibrium points or periodic orbits, and their stability, it is more convenient to have an autonomous system, or at least periodic. In this case, we can transform the system to an autonomous one through the change sin(θ) = tanh(s). After straightforward computations, with the new variables (θ, z, w = z ), the equation (2) becomes the following autonomous system (see [1] for details) where = d ds denotes the derivative with respect to s, and The new variable θ can be interpreted as the time: on one hand, observe that R = 1/ cos 2 θ, therefore, each value of θ provides also the distance between de primaries. On the other hand, when s ∈ (−∞, +∞), θ varies in (−π/2, π/2), and the limits s → ±∞ correspond to θ → ±π/2. Notice that θ > 0, so it is an increasing function: when θ is close to −π/2, the primaries are far apart, then they go close, with the minimum distance when θ = 0, and separate as they escape through their parabolic orbits when θ tends to π/2. We also remark that, if there exists any equilibrium point (as we will see in the next section), it must satisfy θ = ±π/2. Mathematically, the system (4) is well defined when θ = ±π/2 (that is when the primaries are "at infinity"), and the corresponding equations are invariant subsets of the system. We call them upper and lower boundary problem respectively. The key point to consider these boundary problems is that although the study of the dynamics "at infinity" has no physical meaning, it will help to understand some behaviours when the two primaries go close, as we will show. More precisely, a natural way to start studying a nonlinear system of ODE is by the computation of the simplest solutions, the equilibrium points.
We will compute them (at θ = ±π/2) an we will also study their stability, which will give rise to the computation of the associated invariant manifolds that do exist for the whole system 4, also when the primaries are separated by a finite distance. The dynamics of these manifolds in the finite regime, will be relevant for the purpose of this paper: the explanation of formation of bridges and tails.
The extended phase space of system (4) is given by We will call the system (4) the global system, and we denote as configuration space the projection of D on to the (x, y) plane.
Notice that the equations (4) have two singularities at the (fixed) location of the primaries z = z 1 and z = z 2 .
The symmetry implies that for any solution of (4) there exists another one symmetric with respect to y = 0 in the configuration plane reversing time.

Equilibrium points, homothetic solutions and the Jacobi function
All equilibrium points of the global system are points belonging to upper and lower boundary problems, that is, all of them correspond to θ = ±π/2. As was shown in [1], the equilibrium points of the parabolic problem in each boundary coincide with the classical five equilibrium points of the circular restricted three-body problem (see [17], for example): three collinear and two triangular. We denote by L + i and L − i , i = 1, ..., 5, the equilibrium points for θ = π/2, and θ = −π/2 respectively. Along the paper, and in order to follow the same notation as in [1] and [4], we label the collinear equilibrium points increasingly with respect their location on the x axis: Figure 3. The study of the stability character of the equilibrium points is achieved through the study of the eigenvalues of the differential matrix of the vector field of (4) at the equilibrium points (see, for example, [11] for a detailed discussion about the study of the linearisation of the equations around equilibrium points). In the parabolic problem, all the equilibrium points have associated positive and negative real eigenvalues. This means that there exist unstable (for positive eigenvalues) and stable (for negative eigenvalues) invariant manifolds. The definitions of these manifolds for a generic equilibrium point L ξ are the following. Let be z(s) a solution of (4). Then: • The unstable manifold is the set of orbits that tend to the equilibrium point backwards in time: • The stable manifold is the set of orbits that tend to the equilibrium point forwards in time: The invariant manifolds are invariant subsets of the phase space.
In [1], the dimensions of the invariant manifolds of the equilibrium points were studied. In Table 1 we show the dimensions corresponding to the equilibrium points on the upper boundary problem L + i , i = 1, . . . , 5. Using symmetry (7), the dimensions of the invariant manifolds associated to the equilibrium points on the lower boundary problem are obtained. For any collinear equilibrium point L + i , i = 1, 2, 3, the unstable manifold is of dimension one. That means that W u (L + i ) consists in a single orbit. The stable manifold is of dimension four, that is, in the phase space D, W s (L + i ) Table 1: Dimension of the invariant manifolds of the equilibrium points in the global system (4) consists in a 4-dimensional set of orbits that tend to the equilibrium point. This is an important fact: a manifold of dimension n − 1 (or of codimension 1) in a space of dimension n divides the space in subregions. When considering the invariant manifolds associated to the equilibrium points, we will show that W s (L + i ), i = 1, 2, 3 separates the phase space in subregions where the orbits have different behaviours, which gives the clue to understand the formation of bridges and tails.
Besides the equilibrium points, the simplest solutions of the global system (4) are the five homothetic solutions connecting the equilibrium points Figure 2: Clearly these five homothetic solutions in the sidereal or inertial reference frame are solutions in which the three bodies (the infinitesimal mass and the two primaries) keep the same configuration all the time: either the three bodies lie in a line (collinear configuration) or they lie at the vertices of an equilateral triangle (triangular configuration).
The parabolic problem also admits a function, similar to the Jacobi constant of the circular restricted three-body problem (see [17]), that we call, by similarity, the Jacobi function: Unlike the circular problem, C is not constant, in general, along the solutions of the global system (4), but it has a piecewise monotone behaviour along the solutions, except at the homothetic solutions, where its value is constant.
Given a value of C, let V 0 (C) = {z | 2Ω(z) = C} be the set called zero velocity curves. Their topology is the same as in the circular restricted three body problem (see [17]). Let C i = C(L ± i ), i = 1, . . . , 5, be the value of C at each equilibrium point. The so called Hill's regions, that is, the allowed regions of motion in the configuration space, for C fixed, can be obtained from (11).
We plot in Figure  We remark the fact that the zero velocity curves change with time: when s ≤ 0, the curves shrink and the Hill's region gets bigger, whereas for s ≥ 0 the curves get bigger and the Hill's region decreases. This is a key factor for the description of the final evolutions of the solutions. The location of the equilibrium points (although they do not live at the same level of C) is also shown.

Final evolutions
We want to describe all the possible behaviours of the solutions of system (4) as time tends to infinity. We call them the final evolutions of a trajectory: future evolution, when time tends to +∞ (forwards) and past evolution, when time tends to −∞ (backwards). Notice that, using the symmetry (7), for any kind of future evolution there exists a trajectory exhibiting the same kind of evolution backwards in time.
Therefore, the same kind of behaviours are obtained both for past and future. We use in the present paper some results given in [4] for µ = 1/2.
The results are valid for any value of µ, and the proofs are similar, so we do not repeat them here.
We will describe the final evolutions in general (forwards and backwards).
Due to the geometry of the zero velocity curves and the nature of the Jacobi function, the final evolutions are rather simple, mainly, escape and capture orbits. Essentially, an escape orbit is a trajectory that goes far apart from both primaries; a capture orbit is a trajectory that gets trapped around one primary (spinning around it) or even collides with the primary. More specifically, when s → ∞ (representing indistinctly ±∞) any solution of system (4) can be classified as : • Capture orbit. It is an orbit that the synodical distance to one of the primaries tends to zero as time tends to infinity: where j can be 1 or 2. That means that the massless particle approaches one of the primaries and keeps around it as time increases. Eventually, the particle can also collide with the primary (recall the singularities of the equations).
All of the capture orbits tend to a collision with one of the primaries (in general at infinity time). Analogously to the definitions of the invariant manifolds associated to an equilibrium point, the set of all trajectories ending in a collision with the corresponding primary forwards/backwards in time will be denoted by W s (m j ) or W u (m j ), j = 1, 2, respectively. They are called the stable/unstable manifolds associated to collision manifold.
• Escape orbit. It is an orbit that keeps away from any primary as time tends to infinity: for a fixed value of δ. Using (3), this means that the inertial distance |Z(t) − Z j (t)| → ∞ as t → ∞ for j = 1, 2, so the particle escapes from any neighbourhood of any primary.
Notice that, regarding future evolutions, the orbits belonging to . . , 5, satisfy the escape condition. By definition, an orbit on the stable manifold of any equilibrium point L + i will tend to that point as s → +∞. In particular, it does not approach any primary. We want to differentiate these orbits from the orbits that escape without any particular configuration. Focusing on future evolutions: -an orbit escapes in collinear configuration with the primaries if it belongs to W s (L + i ), i = 1, 2, 3; -an orbit escapes in equilateral triangular configuration with the primaries if it belongs to W s (L + i ), i = 4, 5; -an orbit escapes with no configuration, otherwise.
The same definitions can be done for past evolutions and W u ( Any trajectory can be classified depending on their past and future evolution (which, of course, may be different). The evolution of the Hill's regions (again, forwards and backwards in time) allows us to obtain a criterion to classify the Table 2: Position of equilibrium points L ± 2 and the value of orbits depending on their final evolution. We state the criterion for future evolutions, clearly the similar criterion can be obtained for past evolutions. For all values of µ, C 1 < C 2 . Therefore, for future evolutions, it is sufficient to monitor if C(γ(s)) > C 2 to be able to classify the orbit. In Table 2 we show the value of C 2 for the values of µ for which we will present some results in Section 3.

Regularization of the equations of motion
In order to study the collision with a primary and to deal numerically with orbits going close to the primaries, we need to remove the singularities r i = 0, i = 1, 2 appearing in the equations of motion (4). In order to do so, we follow McGehee's ideas (see [14] for the elliptic RTBP and references therein) to obtain the so called regularized system of equations.
The regularization removes one singularity, r 1 = 0 or r 2 = 0 at a time. We describe here the procedure to remove the collision with m 1 . We perform the following changes of variables: 1. We move the selected primary m 1 to the origin: with x = X 1 + µ, y = X 2 , x = X 3 and y = X 4 .
3. We now consider with 4. Finally we introduce a new time variable τ by ds dτ = r 3/2 .
After the computations to implement the above changes of variables and time, system (4) becomes where˙= d/dτ and r 2 = √ r 2 + 2r cos δ + 1. We notice that this new system has only the singularity r 2 = 0.
We proceed in a similar way to remove the singularity r 2 = 0 from system (4), so the new associated system has only the singularity r 1 = 0. Along the integration of any trajectory, for a given initial condition, we integrate numerically system (4) unless the particle is in a neighbourhood of one of the primaries, where we apply the changes of variables and time and integrate the corresponding regularized system of equations.

Dynamics of the parabolic problem
The dynamics of the parabolic problem can be understood focusing on its final and past evolutions and looking for the heteroclinic connections, that is, orbits that connect two invariant objects, such that, two equilibrium points (as the homothetic solutions) or a primary and an equilibrium point (therefore, an orbit belonging to W u (m j ) and W s (L + i )). See [4] for details.
The existence of few different types of final evolutions becomes enriched due to existence of the homothetic solutions. Since they belong to W u (L − i ) ∩ W s (L + i ), they are a natural way to transport the dynamics near the lower boundary problem to the upper one (see Figure 2). Moreover, in the case of the collinear equilibrium points, the invariant manifolds of co-dimension 1 (the orbits that escape, forwards or backwards in time, with linear configuration) behave as a frontier and divide the phase space in regions of capture or escape with no configuration. We will see in Section 3 that these frontiers are also responsible for the existence of bridges and tails.
In order to emphasize the importance of the invariant manifolds of codimension 1 associated to the collinear equilibrium points, and how they separate the different type of orbits, we will reproduce some results from Section 4.2 in [4] (µ = 0.5 in that paper) for values of µ = 0.5. More precisely, we take initial conditions at θ = 0 in the plane (x, y ), that is, with y = x = 0.
These orbits are symmetric with respect to θ = 0, so they have the same final evolution forwards and backwards in time. We integrate these initial conditions forward in time and classify the orbits using the C-criterion. Notice that due to the fact that some orbits can have a close encounter to a primary, a binary collision regularization is performed in order to continue the integration.
In Figure 4 the detailed regions of escape with no configuration (in white) and capture around m 1 (in red) or around m 2 (in blue), on the (x, y ) plane, are shown for different values of µ.
Clearly, according to the possible final evolutions described in Subsection 2.3, the frontiers between escape with no configuration and capture regions must correspond to orbits that escape with a certain configuration, that is, they belong to W s (L + i ) for a certain equilibrium point. As mentioned, the only stable manifolds of co-dimension 1 are those of L + i , i = 1, 2, 3. The frontier between two regions of capture orbits around different primaries correspond to W s (L + 2 ), between an escape with no configuration region and a capture around m 1 is W s (L + 3 ) and between an escape region with no configuration and a capture around m 2 is W s (L + 1 ). To illustrate that, we consider, for µ = 0.3, the region in the (x, y ) orbit a escapes with no configuration, whereas orbit b is captured by m 1 . Therefore, by continuity there must be one orbit tending to L + 3 . That is, the frontier between white and red regions must be orbits belonging to W s (L + 3 ). Analogously, orbits c and d are captured by m 1 whereas orbit e is captured by m 2 . Again by continuity there exist an orbit that tend to L + 2 . Thus, the frontier between red and blue regions must be orbits belonging to W s (L + 2 ). The same reasoning can be done with orbits f (capture m 2 ) and g (escape with no configuration): in between the frontier corresponds to W s (L + 1 ). From now on, the orbits that escape with no configurations will be simply called escape orbits.

Generating bridges and tails
Our goal in this section is to show that the stable invariant manifolds associated to the collinear equilibrium points, W s (L + i ), i = 1, 2, 3, and the unstable invariant manifold associated to the collision with one primary, W u (m i ), i = 1, 2, are responsible for the existence of bridges and tails. We perform a general and broad exploration considering big sets of initial conditions around a primary and we classify them depending on the final evolution forwards in time after the close encounter of the primaries using the C-criterion. We show and comment the results obtained for µ = 0.5, and then for µ = 0.5.
The general exploration is performed as follows. First, we fix a negative value of θ = θ 0 (a time before the close encounter) and a value of C ≥ C 2 .
For this value of C we know that the Hill's region has a bounded component around m j , j = 1, 2 (see Figure 3). Next, we fix a primary, m j and take a circle centred at that primary and radius r c such that it is contained in the bounded component of the Hill's region. Then, we generate a set of initial conditions (x, y, x , y ) around the primary m j as: where x m1 = µ and x m2 = µ − 1, α, β ∈ [0, 2π] and v is obtained from (11). For any given α and β, we take the corresponding initial condition and we follow its trajectory, forward in time. Applying the C-criterion each trajectory will be classified as a collision orbit with one of the primaries or an escape orbit.
For each fixed value of C and r c , we take N equally spaced values of α and β. We show the results in two different ways: • Classification plots. We plot in the (α, β) plane each point in different colours depending on the final evolution: red -captured by m 1 , bluecaptured by m 2 or white (blank) -escape.
• Snapshots. We plot the location of the particles in the configuration inertial frame (X, Y ) for different values of time θ > 0. The code of colours red and blue is the same one. The orbits that escape are plotted in black in these figures.
We will see how, considering different regions in the classification plots (and also for different values of r c and C), the snapshots show the existence of bridges and tails.
Several comments regarding such initial conditions and the classification should be made. We observe that any negative value of θ 0 might be taken. The closer the value of θ 0 to zero, the smaller the distance between the primaries. We do not want to start too far from the closest passage between the primaries, nor too close. Moreover, any particle at a distance r c of the primary in synodical coordinates, is at a distance R c = r c R of that primary in the inertial frame, where R = 1/ cos 2 θ 0 is the inertial distance between the primaries. Thus r c represents the ratio between the distance of the particle to the primary and the distance between the two primaries. We consider θ 0 = −π/4, so R = 2, and r c ∈ (0, 0.2]. Then, we have a set of initial conditions that forms an annulus that spreads around the primary up to 20% the distance between the two primaries in the inertial frame. We also remark that the strategy to take such initial conditions inside the bounded component of neighbourhood the Hill's region guarantees, using the Ccriterion for past evolutions, that all the trajectories considered tend to collision with m j backwards in time. That is, all the trajectories considered belong to the unstable manifold of the collision with m j , W u (m j ).

Results for equal masses
We consider the case µ = 0.5, that corresponds to µ = 0.5. We take a set of test particles around one primary as in (19). Since both primaries have the same mass, it is enough to do the exploration only for one primary. We consider in this section all the particles leaving a neighbourhood of m 1 . We now focus on the frontier separating the different coloured regions. In Section 2.6, we showed that the stable manifolds of L + i , i = 1, 2, 3 separate different types of orbits and the points on the boundary between different coloured regions (see Figure 4 and 5) precisely belong to W s (L + i ), for a suitable i = 1, 2, 3. In a similar way now, we have that the points on the frontier separating the different coloured regions that belong to W u (m 1 ), also belong to: • W s (L + 1 ) if the boundary curve separates a escape region (white) from a region of collision to m 2 (blue); • W s (L + 3 ) if the boundary curve separates a escape region (white) from a region of collision to m 1 (red); • W s (L + 2 ) if the boundary curve separates two different regions of collision (blue and red). Therefore, the curves that separate two different regions belong to W s (L + i )∩ W u (m 1 ), that is, they are heteroclinic orbits connecting a collision with m 1 and an equilibrium point.
Next we show that these heteroclinic orbits provide simple mechanisms to explain tails and bridges. To do so, let us focus on selected ranges in the (α, β) plane where two different coloured regions appear. We choose two zones: a red region (collision to m 1 ) surrounding a white one (escape region) (see the box in Figure 7, left), and a blue region (collision to m 2 ) surrounded by a red one (collision to m 1 ) (see the box in Figure 8, left). We also compute, using bisection method, the points on the frontier (thick black points). Regarding the first zone, we focus on the separation points between escape and capture with m 1 -the black curve in Figure 7, left -, which correspond (backwards and forwards in time) to heteroclinic orbits between m 1 and L + 3 . Close to these heteroclinic orbits, we have escape orbits and orbits such that the particle remains captured by m 1 . This is the mechanism that provokes the apparition of a tail. In Figure 7, right, we plot a snapshot at θ = 1.2 of the location of the particles ((X, Y ) in the inertial frame) with initial conditions (α, β) at the selected zone (the close encounter between the primaries has already taken place at θ = 0). We observe the set of particles -red and black colourthat spread on a large region of the (X, Y ) plane with the shape of a tail. As time tends to infinity (θ → π 2 ), some of them must remain around m 1 , others must escape, and others must tend to the L + 3 point (the thick curve in the right plot).
Concerning the second zone, the points belonging to the frontier between the two capture regions -blue and red (the thick curve in Figure 8, left)-, correspond (integrating backwards and forwards in time) to heteroclinic orbits between m 1 and L + 2 . Similar to the previous case, as close to these heteroclinic orbits as desired, we have orbits such that the particle is captured by m 2 and orbits such that the particle is captured by m 1 , when increasing the time (or θ). This is the mechanism that provokes the apparition of a bridge. In Figure 8, right, we plot a snapshot at θ = 1.2 of the location of the particles ((X, Y ) in the inertial frame) with initial conditions (α, β) at the selected region. We observe the set of particles connecting both primaries with the shape of a bridge. As before, when time tends to infinity (θ → π 2 ), some of them must remain around m 1 , others must transfer to a neighbourhood of m 2 and remain around it, and others must tend to the L + 2 point (the thick curve in the right plot is close to L 2 for θ = 1.2).
To show a bridge and a tail all together, we choose a suitable zone of initial conditions, for example (α, β) ∈ [5, 5.5] × [3,4], r c = 0.00405 ( 1/100 times the maximum r c allowed by the zero velocity curve), and C = 8. In Figure 9 we plot the zone in the (α, β) plane considered (top left) and three different snapshots of the positions (X, Y ) of each particle at θ = 0, θ = π/8 and θ = 1.2. As explained above, the different coloured layers of initial conditions will give rise to bridges and tails. We just remark that, due to the red thin layers between the white and blue regions in the classification plot shown in Figure 9, top left, the tail will be formed behind m 1 , although in Figure 9, bottom right, it is apparently formed around m 2 . As time increases, for θ > 1.

Results for different masses
Although the model considered is academic, we want to consider the case of galaxies of different mass, that corresponds to consider the primaries with different masses, that is µ < 0.5. We do not intend to take particular values of the mass parameter, associated to real pairs of interacting galaxies, but simply we want to check that the mechanism that explains the apparition of tails and bridges also applies for any value of µ. We have taken a set of test particles around both primaries as in (19) for different values of µ and we have encoun- tered heteroclinic orbits connecting the primaries m j , j = 1, 2 with the collinear equilibrium points. The distinction with the case µ = 0.5 is that the primaries have different masses so we have to consider the influence of each primary separately and compute the corresponding proportions of orbits that tend to capture to each primary or escape.
Just as an illustration of the simulations done, in Figure 11 we show a snapshot for µ = 0.3 where a bridge is clearly apparent, and a snapshot for µ = 0.1 to show a bridge and a tail.
Following the procedure to obtain Figure    the initial conditions for θ = −π/4) -we will say that the particles leave m 1 or m 2 , and of course, backward in time, they collide with m 1 or m 2 -. The solid and dashed lines correspond to orbits that leave m 1 and m 2 respectively. The colours orange, violet and green correspond to µ = 0.5, µ = 0.3 and µ = 0.1 respectively (labels in the plots have been added for black-white print).
As we can see from the left plot, for C close to (and bigger than) C(L 2 ), the smaller the value of µ, the bigger the number of orbits captured by m 1 (assuming that the particles leave m 1 or m 2 ). However, when the particles leave m 2 , the tendency in the proportion of capture orbits by m 1 is inverted when C increases, which seems reasonable because with bigger values of C, the Hill's region around m 2 shrinks with decreasing µ and there is a less possibility to leave the neighbourhood of m 2 .
From the centre plot, an opposite behaviour can be observed for the proportion of capture orbits by m 2 . For C close to (and bigger than) C(L 2 ), this proportion decreases with µ (when leaving m 1 or m 2 ), but the tendency is inverted when C increases, that is the proportion increases with µ, just taking into account the particles that leave m 2 .
Finally, taking into account both proportions of capture by the primaries, we obtain the right plot for the escape orbits. We observe, in particular, that given µ and for (suitable) big values of C, there are no escape orbits.

Discussion
At this point, it seems appropriate to discuss how our results compare with previous ones in the literature.
As mentioned in the Introduction, several papers analyse the observed bridges and tails for pairs of interacting galaxies. These papers use different models and a first topic in this Discussion Section is related to the model we have considered, the parabolic RTBP. Two remarks must be done concerning this model: (i) we assume that each galaxy is modelled as a disk of non-interacting particles orbiting around a central mass point (a primary), whereas the two mass points move in parabolic orbits with respect to their centre of mass. In particular trajectories described by any particle that pass very close to a primary or even collide with it may take place. A regularization of the system of ODE must be carried out as an strategy to deal with the numerical integration of the system close to or at a singularity (collision with a primary). (ii) Of course the parabolic RTBP is a very simple model, but we emphasize that our purpose in this paper is to show a mechanism, using dynamical systems tools and more particularly using invariant manifolds, that explains the formation of tails and bridges in this simple model.
As a second topic, and concerning the references mentioned in the Introduction, we will focus on [18] paper, since our interest for the problem starts with their work, where the authors show that bridges and tails appear when the encounter of two galaxies is also modeled by the parabolic problem.
The authors consider two situations: equal masses (which corresponds to µ = 0.5) and when the big primary is four times bigger than the small one (that is µ = 0.2). Once µ is fixed, they take a bunch of particles -located in different rings-around one primary (typically the big one) and far from the other primary, which is bare. In this situation, the dynamics around the primary can be modelled by a two body problem (the influence of the second primary is a pertur- We have seen that a dynamical explanation for bridges and tails comes from the existence of heteroclinic connections between collision with the primaries and the equilibrium points L + i , i = 1, 2, 3. Here we want to reproduce some explorations of Toomre and Toomre's paper in order to show that the initial conditions that they considered are close to these heteroclinic connections, so this is the reason why they see bridges and tails.
More precisely, we perform the explorations in two ways: 1. Method 1 (TT 1 ): we recall that in the Section 3, we have obtained the classification plots (Figures 6, 7 and 8), considering all the possible initial conditions in position (varying the synodical radius r c ) and velocity for a given C. In this first simulation, we want to restrict, among this whole set of initial conditions, to those synodical initial conditions (x, y, x , y ) as in (19), at θ = −π/4, such that the sidereal velocity (X , Y ) is perpendicular to the position vector (X, Y ), centred at the big primary. To satisfy this restriction, for each value of α, only two values of β are admissible: where the ± sign corresponds to a direct/retrograde orbit.
This means that the initial conditions are at the apoapsis or periapsis of their orbits around the primary, although we do not ensure that the orbit is circular (thinking in a two-body problem m 1 plus a particle) because the modulus of the synodic velocity v is obtained from the fixed value of  At θ 0 = −π/4, we take two sets of sidereal initial conditions in rings of radius R c centred at the primary, that is, (X, Y, X , Y ) as: such that V 2 = (1 − µ)/R c for j = 1 or , that is, velocity to ensure initial conditions on a circular orbit and β = α±π/2 (so the position with respect the primary (X mj , Y mj ) and the velocity vectors are perpendicular); as before, the ± sign refers to a direct/retrograde orbit, -called from now on direct and retrograde cases respectively-.
For each sidereal initial condition, we compute the corresponding synodic initial condition (x, y, x , y ) that can be written as in (19) with r c = R c /2, although now all of them have a different value of C Concerning the results obtained taking rings of particles around the big primary m 1 , in Figure 14 we show the classification plots for the direct case. We see that only collision orbits with m 1 and escape orbits appear, and that as µ decreases, all the initial conditions correspond to collision orbits.
Recall that the smaller the value of µ, the bigger the mass of m 1 . Therefore, direct orbits only contribute to the formation of tails if µ is not too small.
Regarding the retrograde case, we observe from Figure 15  with them concluding that retrograde initial particles provide both tails and bridges.
• Taking into account different mass galaxies, and initial rings of particles around m 1 , Toomre and Toomre just consider two cases: µ = 0.5 and µ = 0.2, but in the latter they only examine the retrograde case. In this paper, we have taken several different values of µ ∈ (0, 0.5], and also both the direct/retrograde cases. We can conclude from the plots in Figures 14 and 15   And we conclude that there are good (clearly visible) bridges even for equal masses.
Concerning the results obtained taking initial rings of particles around the small primary m 2 , in Figures 16 and 17 we show the classification plots for direct and retrograde orbits, respectively. In the case µ = 0.5 we obtain the symmetric plot of Figure 14 and 15 top left, due to the symmetry of the problem. As µ decreases, more diversity appears in both direct and retrograde orbits, and in particular the set of orbits that tend to m 1 grows. This can be explained simply by the fact that as µ decreases, m 1 gets bigger. It is worth mentioning that in Toomre and Toomre's simulations, they assert that for µ = 0.2 (retrograde case) -the only case they consider-, there do not appear enduring bridges. However, as we have seen from our simulations, for µ ≤ 0.3, in the direct case and for all µ ∈ (0, 0.5] in the retrograde one, there do appear bridges, since coloured red-blue regions are obtained.

See Figures 16 and 17.
A final comment is that both in Toomre and Toomre's paper and in method 2 of our computations, only initial rings assumed to describe circular orbits are considered. Of course, other effects as eccentricity and inclination might be taken into account (actually they consider inclination and we do not). For realistic models, of course, such effects are relevant, but we insist that in this paper, the purpose was to present a mechanism that explains the formation of bridges and tails, using dynamical systems tools and a very simple model.
The next natural step is to take other more realistic models and compare and contrast the new simulations obtained (concerning bridges and tails) with the morphology of particular observed real interacting galaxies. But this will be the goal of a future paper.

Conclusions
Although some observations are published showing the existence of bridges and tails resulting from interacting galaxies, as far as we know, this is the first attempt to explain their formation applying invariant manifolds.
In this paper, we have considered a very simple model, the planar parabolic RTBP, where each galaxy is modelled as a disk of non-interacting particles orbiting around a central mass point (a primary) whereas the two primaries move in parabolic orbits with respect to their centre of mass. A regularization of the system of ODE has been carried out in order to numerically integrate this system allowing close passages or even collisions between a particle and the central mass.
A dynamical mechanism that explains the formation of bridges and tails after a close approach of the two galaxies is analysed. Such an explanation is due to the existence of heteroclinic orbits between the collinear equilibrium points and the collision manifold associated to the primaries. The fact that such invariant manifolds separate different final dynamical evolutions is directly responsible for the formation of tails and bridges. We have done some numerical simulations where the tails and bridges are clearly seen. In addition, we have classified whole sets of initial conditions for particles (located in rings) that were orbiting one primary, in a direct or retrograde direction, and after a close encounter between the two galaxies, may remain captured by the same primary, or jump and remain captured by the other one or escape.
Roughly speaking, we can conclude that, on one hand, tails are formed either from particles that were orbiting the small primary or from those orbiting the big one if the mass parameter µ is big. On the other hand, bridges are formed either from particles that were orbiting the small primary of from those orbiting the big one in a retrograde direction for big values of µ.
Finally, the mechanism described has been applied to explain previous simu-