Revealing roaming on the double Morse potential energy surface with Lagrangian descriptors

In this paper, we analyse the phase space structure of the roaming dynamics in a 2 degree of freedom potential energy surface consisting of two identical planar Morse potentials separated by a distance. This potential energy surface was previously studied in Carpenter B K et al (2018 Regul. Chaotic Dyn. 23 60–79), and it has two potential wells surrounded by an unbounded flat region containing no critical points. We study the phase space mechanism for the transference between the wells using the method of Lagrangian descriptors.


Introduction
The purpose of this paper is to reveal the phase space mechanism for roaming in the double Morse potential energy surface using the method of Lagrangian descriptors. The roaming mechanism for chemical reactions was first studied in this model in [1].
The roaming mechanism for chemical reactions was reported in [2]. This paper inspired and motivated further work on the topic that is continuing. Some recent reviews of the topic are [3][4][5][6][7].
Following [1], a concise description of the roaming mechanism that seems to capture all such reactive phenomena attributed to roaming to date is the following. An energised molecule, which appears on the verge of dissociating into two fragments, instead exhibits behaviour 1 Author to whom any correspondence should be addressed.
Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
where one fragment begins to rotate about the other in the flat region of the potential energy surface, culminating in the two fragments re-encountering each other in a 'reactive event'. This reactive event may be constitutionally identical to one for which a familiar transition state associated with an index one saddle can be identified, but the roaming reaction path is observed to avoid the intrinsic reaction path.
The double Morse potential energy surface allows us to analyse key feature of the reaction dynamics which appear to underlay the roaming mechanism. In particular, the double Morse has a flat 'roaming region' where it can be rigorously proven that there are no saddle points. For explaining reaction mechanisms in the context of the potential energy surface index one saddles have been a fundamental building block. This has motivated efforts to locate index one saddle points that can be attributed as the mechanism for roaming ('roaming saddles').
However, the work in [1] showed that roaming mechanism in the double Morse potential occurs without the influence of saddle points. In particular, it was demonstrated that the roaming was governed by phase space structures governing dynamics in the flat region of the potential energy surface. These phase space structures are normally hyperbolic invariant manifolds (NHIMs) [8][9][10][11][12], which in the case of 2 degree of freedom Hamiltonian systems are unstable periodic orbits [13].
The main issue of interest in [1] was the motion of trajectories between the two wells and the intermediary dynamics in the flat region. It was conjectured that three types of unstable periodic orbits mediated the dynamics. The description of roaming given in [1] was based on the behaviour of trajectories initiated on periodic orbit dividing surfaces associated with the three types of periodic orbits. The trajectory behaviour confirmed the conjecture. However, this work did not examine the phase space structures that determine the different types of trajectory behaviour. Such an analysis would involve an understanding of the geometry of the stable and unstable manifolds of the periodic orbits in phase space, similar in spirit to the analysis that has been carried out in [14,15].
In this paper, we will determine the phase space pathways for the roaming mechanism in the double Morse oscillator using the method of Lagrangian descriptors. Lagrangian descriptors are a trajectory diagnostic that has proven useful for revealing phase space structures in dynamical systems. The method was originally developed to reveal fluid flow structures in Lagrangian transport problems (hence the name Lagrangian descriptors) [16]. However, in recent years, the ease of implementation and interpretation of the method has applied in fields far beyond fluid mechanics. In particular, in theoretical chemistry the method has seen a diverse collection of applications to problems in chemical reaction dynamics in recent years, see, e.g., [17][18][19][20][21][22][23][24][25].
The method facilitates a high-resolution approach for exploring high dimensional phase space with low dimensional sets. It can also be applied to both Hamiltonian and non-Hamiltonian systems [26] as well as to systems with arbitrary, even stochastic, timedependence [27]. Moreover, Lagrangian descriptors can be applied directly to trajectory data sets, without the need for an explicit dynamical system [28].
The 'classical' Lagrangian descriptor field is computed by selecting a region of phase space and then choosing a grid of initial conditions for trajectories in this set. Points in the set are assigned a value according to the arc length of the trajectory starting at that initial condition after integration for a fixed, finite time, both backward and forward in time (all initial conditions in the grid are integrated for the same time). The idea underlying Lagrangian descriptors is that the influence of phase space structures on trajectories will result in differences in arc length of nearby trajectories. This has been verified and quantified in terms of the notion of singular structures in the Lagrangian descriptor fields, which are easy to recognise visually [26,29].
Trajectories are the elementary objects that are used to explore phase space. In fact, phase space structure is built from trajectories and their behaviour. For high dimensional phase space this approach is problematic and prone to issues of interpretation since a tightly grouped set of initial conditions may result in trajectories that become far with respect to each other in high dimensional phase space. The method of Lagrangian descriptors takes a new and different approach by encoding the phase space structure in the initial conditions of trajectories, rather than the precise location of their futures and pasts (after a specified amount of time). Hence, a low dimensional set of phase space can be selected and sampled with a grid of initial conditions of high resolution. Since the phase space structure is encoded in the initial conditions of the trajectories, no resolution in the chosen slice of phase space is lost as the trajectories evolve in time.
This paper is outlined as follows. In section 2, we recall the main features of the double Morse potential energy surface and resulting Hamiltonian system as described in [1] that are relevant to our analysis. In section 3, we define the features of the double Morse potential energy surface that are relevant to our analysis. In section 4 we describe Lagrangian descriptors. In section 5, we show how Lagrangian descriptors can be used to compute periodic orbits and their stable and unstable manifolds. In section 6, we describe the phase space pathways for roaming using the stable and unstable manifolds of the periodic orbits.

The double Morse potential energy surface
The 2 degree of freedom Hamiltonian system that we consider is the sum of two identical planar Morse potentials [30] separated by a distance 2b: and kinetic energy of the form: Therefore the Hamiltonian is: Following [1], the parameters for the Hamiltonian are chosen as follows: The parameter D e is the depth of the Morse potential well. Therefore for total energies above D e trajectories can dissociate, i.e. they become unbounded, and the value of this parameter is referred to as the dissociation threshold. Taking the units of energy to be kcalmol −1 , and units of distance to be Å, the parameters of the single Morse potential functions in equation (1) correspond, roughly, to the vibration of a CH bond. Chemically, this potential could be considered  to represent an atom interacting with a diatomic molecule, or, mechanically, a rigid dumbbell interacting with a particle. The equilibrium bond length between the two atoms of the diatomic is set by b, and the bond length of the third atom to one of the atoms of the diatomic is set by r e .

Features of the double Morse potential energy surface
The figures 1 and 2 show plots of the potential energy in colour scale and equipotential lines. At the origin, (x = 0, y = 0) the potential energy has an index one saddle. The potential energy evaluated in the saddle point is V(0, 0) = E s . This index one saddle point exists for all values of b, and E s is below D e in this example, see figure 3. The double Morse potential energy function goes exponentially fast to zero in the asymptotic flat region. This almost flat region plays an essential role in the roaming phenomenon, and it is called the roaming region.

Lagrangian descriptors and 2 degree of freedom Hamiltonians systems
The 'classical' Lagrangian descriptor is based on the arc length of a trajectory. In the present article, we are going to use a variant based on the integral of the square of the infinitesimal arc length proposed in [26]. Moreover, the present variant considers two modification in the definition of the Lagrangian descriptor. The purpose of the first modification is to differentiate between the phase space structures defined with the propagation of the trajectories forward and backwards in time. The second modification in the definition of Lagrangian descriptor allows the different time propagation of the trajectories for different initial conditions in order to study the phase space structures contained in a specific region of the phase space. We consider these modifications in order to improve the visibility and reveal the important phase space structures in the present study of roaming and dissociation. The definition of this variant of the Lagrangian descriptor is as follows.
Let us consider the autonomous ordinary differential equations where v(x, t) ∈ C r (r 1) in x and continuous in time. The definition of Lagrangian descriptor depends on the initial condition x 0 = x(t 0 ), on the time interval [t 0 + τ − , t 0 + τ + ], and takes the form, where the over dot symbol represents the derivative with respect to the time t and τ + 0 and τ − 0 are freely chosen parameters. In the 'classical' Lagrangian descriptor the magnitude of τ + and τ − are equal for all points x 0 . In the present variant of Lagrangian descriptor, the values of τ + and τ − can change between different initial conditions. This modification allows us to stop the integration once a trajectory leaves a specific region in the phase space and reveal only the structures related contained in this region. The phase space of a 2 degree of freedom Hamiltonian system has 4 dimensions. Considering the conservation of the energy, it is possible to represent the dynamics of the system in the three-dimensional constant energy manifold. In the constant energy manifold, it is possible to visualise the dynamics and identify the essential structures to understand the dynamics.
The periodic orbits are essential structures to understand the dynamics of the system in the constant energy manifold. Around a stable periodic orbit, the KAM-tori confine the trajectories in the region defined by the tori. In contrast, the dynamics in a neighbourhood of an unstable hyperbolic periodic orbit has a different nature. Its stable and unstable manifolds intersect in the hyperbolic periodic orbit and direct the trajectories in the neighbourhood. The definition of the stable and unstable manifolds W s/u (Γ) of the periodic unstable orbit Γ is the following, In other words, the stable manifold W s (Γ) is the set of trajectories that converge to the periodic orbit Γ as the time t goes to ∞. The definition for the unstable manifold W u (Γ) is analogous, that is, it is the set of trajectories converging to the periodic orbit as the time t goes to −∞.
A 2 degree of freedom Hamiltonian system, W s/u (Γ) has 2 dimensions and form impenetrable barriers in the constant energy manifold. The segments of the stable and unstable manifolds direct the transport between different regions [31].
Another important property of the stable and unstable manifolds related to the chaotic dynamics is that, if a stable manifold and an unstable manifold intersect transversally at one point, then an infinite number of transversal intersections between them exist. The structure generated by the union of the stable and unstable manifolds is called tangle and defines a set of tubes that direct the dynamics in the phase space. The trajectories in a tube never cross the boundaries of a tube. This fact is a consequence of the uniqueness of the solution of the ordinary differential equations.
The Lagrangian descriptors are convenient tools to explore the phase space structure, in particular, to find stable and unstable manifolds of periodic orbits. To understand the basic idea that underpins the detection, let us consider the behaviour of the trajectories in a neighbourhood of stable manifold W s (Γ). The trajectories in W s (Γ) converge to the periodic orbit Γ, and the nearby trajectories in the neighbourhood of W s (Γ) have similar behaviour for a finite interval of time. After this interval of time, the trajectories move away from the stable hyperbolic periodic orbit Γ following the unstable manifold W u (Γ). This different behaviour generates the singularities in the Lagrangian descriptors [26].
The variant of Lagrangian descriptor in the equation (6) has negative sign between the two integrals. This sign allows differentiating between the stable and unstable manifolds on the Lagrangian descriptor plots. The set of singularities on integral with τ + shows the stable manifolds and the set of singularities on integral with τ − the unstable manifolds. The superposition of the two patterns of singularities reveals the tangle between the stable and unstable manifolds.

Periodic orbits and their stable and unstable manifolds
In order to show the dynamics of this 2 degree of freedom Hamiltonian system in the regime of dissociation and roaming, let us consider a fixed value of the total energy above the dissociation energy threshold, E = 101 > D e . In this case, the trajectories close to the wells can cross the interaction region and escape to infinity. The trajectories that start in the interaction region are sensitive to initial conditions. That means that two trajectories that start arbitrarily close can finish in a different direction in the asymptotic region. This phenomenon is called chaotic scattering and is a kind of transient chaos [32,33].
The chaotic scattering is a common phenomenon in the open Hamiltonian system. Detailed studies of chaotic scattering in Hamiltonian systems with 2 and 3 degrees of freedom are in the recent literature [34][35][36][37]. The tangle between the stable and unstable manifolds generates rich dynamics in the interaction region. However, the dynamics is simple in the asymptotic region, where the trajectories converge to straight lines.
The double Morse potential energy V satisfies the asymptotic condition Vr 2 → 0 when r → ∞, which implies that the linear momenta of the trajectories that reach the asymptotic region converge to a constant value [38].
The trajectories associated with dissociation reactions start in one well and finish in the asymptotic region. The trajectories associated with roaming reactions start in one well, travel outside the well in the roaming region where the potential energy is almost flat, and reach the other potential well.
In order to define the transport between the two potential wells and the asymptotic region, the configuration space is divided into three regions, see figure 4. Two circumferences around each potential well define the regions A and B, and another external circumference around the two wells defines the asymptotic region C where the trajectories are close to straight lines and escape to infinity.
A relevant conclusion in the work [1], is that the roaming and dissociation trajectories in the system are related to three types of periodic orbits. Let be Γ 1 , Γ 2 , and Γ 3 periodic orbits representatives of each type. The projections of the three periodic orbits in configuration space is shown in figure 5. The projection of periodic orbit Γ 1 encircles the interaction region and  the two potential wells. The projection of Γ 2 also encircles the two wells but self intersects at the origin. The projection of Γ 3 encircles just one well.
In the configuration space, the projection of periodic orbits Γ 2 , Γ 3 are close, but it is clear that the orbits are separated in the constant energy manifold, see figure 6. The three orbits are very close to each other in a neighbourhood of the point (x = 8.45, y = 0, p x = 0).
As a consequence of the symmetries x → −x and y → −y of the potential energy V(x, y), there are 4 orbits of type I, 2 orbits of type II, and 4 orbits of type III in the phase space. It is necessary to consider only one element of every type for the numerical calculations.
A common method to find the periodic orbits is based on the Poincaré map defined on a two-dimensional Poincaré surface. The Poincaré map is a function from the Poincaré surface to itself obtained by following trajectories from one intersection of the Poincaré surface to the next. The fixed points of the Poincaré map are the intersections between the periodic orbits and the Poincaré surface. In this system, a natural choice of Poincaré surface is the canonical plane x-p x with y = 0 due to the symmetry y → −y of the potential energy. The three periodic orbits Γ 1 , Γ 2 , and Γ 3 are hyperbolic and considerably unstable. In order to calculate the eigenvalues and eigenvectors of the linear approximation to the Poincaré map around the fixed points associated with the three types of periodic orbits is necessary to use multi-precision integrators. These calculations are done with the Taylor integrator method implemented in Julia [39]. The periods of the periodic orbits and their associated eigenvalues are in table 1. The product between the eigenvalues is close to 1, so the calculations agree with the Hamiltonian properties of the system. It is possible to construct the stable and unstable manifolds of the periodic orbits using their corresponding eigenvalues and eigenvectors. However, the construction and visualisation of the stable and unstable manifolds are not simple tasks because of the size of the eigenvalues, λ k,max 1 and λ k,min ∼ = 0 for k = 1, 2, 3, and because the Poincaré map requires that the trajectories intersect the Poincaré section again. For example, the trajectories in the unstable manifold of the periodic orbit Γ 1 that go directly to the asymptotic region intersect the Poincaré section only for a finite number of iterations before to converging to straight lines. In general, it is not possible to obtain a global visualisation of the tangles between the stable and unstable manifolds using the Poincaré map.
An appropriate approach that allows an easier visualisation of the stable and unstable manifolds even in this situation is the use of indicators for chaotic regions, like fast Lyapunov indicators (FLI) [40], mean exponential growth factor of nearby orbits (MEGNO) [41], the smaller (SALI) and the generalized (GALI) alignment indices [42], delay time functions, scattering functions [37], and Lagrangian descriptors. The Lagrangian descriptors contain information about the stable and unstable manifolds that intersect the set of initial conditions and do not require that the trajectories intersect any particular predefined set again, as mentioned in section 4. It is important to notice that the trajectories should be integrated for long enough time to reach a neighbourhood of the geometric structures that generate different abrupt behaviour. In the other case, the chaotic indicator does not show the structure clearly, even if the set of initial conditions intersect the structure.
The calculations of the Lagrangian descriptors are done with the classical ninth order Vernet integrator implemented in Julia [43]. In this system, the Lagrangian descriptors require standard double-precision to reveal segments of the invariant manifolds. The computational cost necessary to see a signature of the stable and unstable manifolds in the Lagrangian descriptor plots is a fraction of the computational cost necessary using the Poincaré map. The reason for this difference is related to the computational cost of the integration with the multi-precision libraries and the fact that the calculation of the Lagrangian descriptor does not require that the trajectories cross the set of initial conditions again.
The figure 5 shows that Γ 1 , Γ 2 , and Γ 3 intersect the x axis in a neighbourhood of the point (x = 8.5, y = 0) and their momentum at the intersection are (p y < 0, p x = 0). In order to find these periodic orbits it is natural to calculate the Lagrangian descriptor evaluated on the canonical conjugate plane p x 0 -x 0 , y 0 = 0, and initial momentum p y 0 < 0. The Lagrangian descriptor evaluated in this set of initial conditions for different times is shown in figure 7. When the size of the interval of integration [τ − , τ + ] is increased, the Lagrangian descriptor reveals intersection of the the stable and unstable manifolds with the set of initial conditions. The intersections Figure 7. Lagrangian descriptor M evaluated on the plane p x 0x 0 with initial y 0 = 0 and initial momentum (p y 0 < 0, p x = 0) for different values of τ + and τ − . To show the symmetry of M and the invariant manifolds for this initial conditions, τ + = −τ − is selected. When the value of τ + is increased, the structure of the tangle between the stable and unstable manifolds is revealed by the abrupt jumps in the value of the Lagrangian descriptor where it converges to a non-differentiable function. The intense blue and yellow curves correspond to segments of the unstable manifolds and stable manifolds respectively. between the stable and unstable manifolds is the invariant chaotic set. The periodic orbits Γ 1 , Γ 2 , and Γ 3 are contained in this invariant chaotic set.
In order to show clearly the peaks generated by invariant stable and unstable manifolds of the three hyperbolic periodic orbits in the Lagrangian descriptor plots and find the periodic  orbits, let us consider the symmetry line y 0 = 0 in intervals that contain the periodic orbits. The Lagrangian descriptor in figure 7 with τ − = −τ + evaluated on this line is identically zero due to the symmetries of the system. For this reason, is convenient to consider the Lagrangian descriptor with τ − = 0 and τ + > 0, the results are in the figure 8. When the integration time τ + is increased enough, the peak generated for the different behaviour of the trajectories reveals the intersections of the stable manifold of the periodic orbit with the x 0 axis. Due to the symmetry of the stable and unstable manifolds with respect to the x 0 axis, this intersection coincides with the intersection with the unstable manifold of the periodic orbit. Finally, the intersection point between the stable and unstable manifolds is the intersection with the periodic orbit with the x 0 axis. It is important to notice that the integration time to generate the peak is different from the period of the orbit. For the periodic orbits Γ 1 and Γ 2 the time necessary to generate the peak is less than their corresponding periods.
Another simple way to find the intersections between the stable and unstable manifolds is to change the negative sign in the definition (6), then M + + M − > 0 and the points where is this quantity is maximal are the intersections between the stable and unstable manifolds.  [35,37].
In figure 10 there are three disconnected regions: two unbounded regions on the extremes and one bounded region in the middle surrounded by a white region. The white region corresponds to points where the momentum p x 0 is not defined for this value of the energy.

Phase space pathways for roaming
An important set in the phase space for analysing the transport in the system is the dividing surface of a periodic orbit [1,44]. The projection of the periodic orbit in the configuration space defines a closed curve and a region inside the curve. The dividing surface is a natural surface to analyse all the trajectories of the system that cross the boundary of this region in the configuration space. To construct the dividing surface in the phase space, let us consider all trajectories that cross the boundary of this region in the configuration space. For every point in the boundary, we consider all the possible momenta compatible with the energy of the system. This procedure generates a set of points in the phase space for every point on the boundary of the region in the configuration space. The union of those sets of points forms a surface in the phase space.
In more detail, the algorithm to construct the dividing surface of the periodic orbit Γ consists of three steps: • Project the periodic orbit Γ on the configuration space.
• Construct the circumference on momentum plane using the equation for every point (x, y) on the projection of the orbit Γ in the configuration space. • Take the union of all these circumferences in the phase space to construct the dividing surface.
The dividing surface constructed in this way has three relevant properties for analysis of the transport: the first property is that the periodic orbit Γ and the orbit with the same projection in configuration space but opposite momentum are contained in the same dividing surface. This fact is a consequence of the time-reversal invariance of the system. The second property is that the periodic orbits used to construct the dividing surface are the boundaries between the regions where the trajectories enter into the phase space region contained by the dividing surface and trajectories that left the same region. The third property, the flux through the dividing surface is minimal. That means if the dividing surface is deformed, the flux through it increases. The demonstrations of these properties of the dividing surfaces of periodic orbits are in [44], and for NHIMs with more dimensions in [45].
In the system, the periodic orbits Γ 1 , Γ 2 , and Γ 3 are not generated by any saddle in the potential energy. Their projections in the configuration space are curves that enclose an area. In the case of unstable periodic orbits associated with saddles points, the projection on the configuration space is a curve that does not enclose any area. The extremes of the projection are the returning points where the periodic orbit has zero momentum. The dividing surface, in this case, is a surface genus 0 (topologically equivalent to a sphere). An example of this kind of periodic orbit and its corresponding dividing surface is in [46].
The projection in configuration space of the orbits Γ 1 , Γ 3 do not have self-intersections; for this reason, the corresponding dividing surfaces are genus 1 (a torus). Another example of this kind of dividing surface is in [47]. The projection of the periodic orbit Γ 2 has an intersection at the origin, the dividing surface, in this case, is a surface genus 2.
A natural parametrisation of the dividing surface is given by the arc length l 0 of the respective periodic orbit in the phase space and the angle between the initial momentum and the x 0 axis, φ 0 = arctan(p y 0 /p x 0 ). The next calculation uses this parametrisation to analyse the trajectories that cross each of the three dividing surfaces.
The fate map is an appropriate tool to classify the trajectories associated with rooming and dissociation in phase space. A fate map shows the origin and destiny of the trajectories, considering the evolution for an interval of time. In this case, the origin and destiny of the trajectories are the regions A, B, and C defined in the previous section. The procedure to calculate the fate map on a dividing surface is the following: consider a point on the dividing surface, integrate forward and backwards the trajectory that crosses this point on the dividing surface until each extreme of the trajectory reaches one of the three regions, and label this point with the corresponding fate. The left side of figure 11 shows the fate map for each dividing surface. Each different colour in the fate map indicates different transport between regions A, B, and C. The fate maps evaluated on the dividing surfaces of the periodic orbits Γ 1 and Γ 2 have all nine possibilities of transport between the regions A, B, and C. In the case for the dividing surface of the periodic orbit Γ 3 , the fate map has only five possibilities.
In the three fate maps there exist red regions (C → A), and light blue regions (A → C). The light blue regions represent the trajectories that dissociate. The green and orange regions,  In the configuration space, all the trajectories that travel from one well to the other well should cross the projection of the periodic orbits Γ 3 and Γ 2 . For this reason, the trajectories associated with roaming cross the dividing surface of Γ 3 and Γ 2 in the phase space.
The trajectories that start in one well and escape to the asymptotic region need to cross the projection of Γ 1 in the configuration space as well. These trajectories are associated with dissociation. The fate map of the dividing surface of Γ 1 shows which trajectories escape to the asymptotic region and which trajectories enter the interaction region.
In order to show the relation between the different trajectories' fates and the stable and unstable manifolds of the unstable periodic orbits Γ 1 , Γ 2 , and Γ 3 we calculate the Lagrangian descriptors on the dividing surfaces for the three hyperbolic periodic orbits. The same segments of the trajectories are used to calculate the Lagrangian descriptors and the fate maps. Like in the calculation of the fate maps, the values of τ − and τ + for each trajectory are the times necessary to reach one of the three regions. Once a trajectory reaches one of the three regions, its integration stops. The results for the Lagrangian descriptors are on the right side of figure 11.
The abrupt jumps in the values of the Lagrangian descriptors are related to the different behaviours of the trajectories in a neighbourhood of the stable and unstable manifolds of the unstable periodic orbits in the phase space. The stable and unstable manifolds are surfaces of dimension 2, and their intersections with a dividing surface are curves. There is a complete match between the boundaries of regions in the fate maps and the abrupt changes in the Lagrangian descriptors plots for each dividing surface. This agreement is a manifestation of the division of the phase space into the roaming and dissociation regions by the stable and unstable manifolds.
There are three different types of regions in the plots of the Lagrangian descriptor evaluated on the dividing surface: blue, yellow, and mixed regions. The colours in the Lagrangian descriptor evaluated in the dividing surfaces are related to the fates. The blue regions in the Lagrangian descriptors plots correspond to trajectories that start in one potential well, cross the dividing surface, and finally escape to infinity. These are trajectories associated with dissociation. The yellow regions correspond to trajectories that start in the asymptotic region, cross the dividing surfaces, and finally reach a potential well. The regions with mixed colours in the Lagrangian descriptor plot without abrupt changes are associated with the trajectories that connect the two wells. The mixed colours are related to other segments of the stable and  unstable manifolds contained in the same set of initial conditions, and more integration time is necessary to visualise them.
The procedure to identify the stable and unstable manifolds associated with the transport from the Lagrangian descriptor and fate map plots is the following. For illustration, consider magnifications of the plots around the point (l 2 = 25, φ 2 = 2.565), close to the periodic orbit Γ 2 , see figure 12. There are fine strips close to Γ 2 on both magnifications. Then consider the vertical line of initial conditions in figures 12(a) and 12(b), and their trajectories used to calculate the fate map and the Lagrangian descriptor. The projection of these orbits in the configuration space is in figure 14. For example, the trajectory on the boundary between the green region (A → B) and blue region (A → C) converges to the invariant manifold W s (Γ 2 ), see the figures 12(a) and 14(a). Continuing with the same procedure for the other boundaries is possible to identify the invariant manifolds associated with the boundaries between regions, see figures 14(b)-(d).

Conclusions and final remarks
The double Morse system is an example of the two centre problem with a rich phase space structure. The periodic orbits type I, II, and III associated with the roaming and the dissociation phenomena are hyperbolic and considerably unstable (λ k,max 1, k = 1, 2, 3) for energies close to the asymptotic value of the potential energy.
Lagrangian descriptors are useful tools for finding the invariant manifolds of considerably unstable periodic orbits in a simple way. The abrupt peaks in the Lagrangian descriptor plots reveals the intersections between the set of initial conditions and the stable and unstable manifolds. The set of trajectories that generate the abrupt peaks in the Lagrangian descriptors plots do not necessarily return to the set of initial conditions. This property of the Lagrangian descriptors allows one to easily study the stable and unstable manifolds in the phase space of open systems. The principle that generates those peaks is based on the different behaviour of the trajectories in the stable and unstable manifolds and the trajectories in their neighbourhoods.
In the case of 2 degree of freedom Hamiltonian systems, it is easy to visualise the structure of the tangle between the stable and unstable manifolds using 2 dimensional sets of initial conditions. The stable and unstable manifolds in the constant energy manifolds have 2 dimensions, then their intersection with the set of initial conditions is a set of curves. In Hamiltonian systems with more than 2 degrees of freedom, the sets of initial conditions necessary to completely visualise the structure of the tangles between stable and unstable manifolds of periodic orbits and NHIMs have more dimensions. Their complete visualisation and analysis is a challenging open problem. Some recent works on visualisation and analysis of phase space structures for 3 degree of freedom Hamiltonian systems are [36,37,[48][49][50][51][52][53].
The dividing surfaces of periodic orbits of type I, II, and III are convenient surfaces to study the transport between the potential wells. The fate maps evaluated on the dividing surfaces of those periodic orbits allow one to analyse the transport between the regions that define the dissociation and roaming. The boundaries of the different regions on the fate map correspond to the stable and unstable manifolds of hyperbolic periodic orbits on the Lagrangian descriptor plots. The stable and unstable manifolds of type I, II, and III of periodic orbits determine the dissociation and roaming in this system.
The stable and unstable manifolds of the hyperbolic periodic orbits type I, II, and III are robust under perturbations due to the large instability of the periodic orbits. Therefore, it is possible to observe a similar mechanism for roaming and dissociation in other systems with two wells in the potential energy when the energy is close to the dissociation threshold.