The phase space geometry underlying roaming reaction dynamics

Recent studies have found an unusual way of dissociation in formaldehyde. It can be characterized by a hydrogen atom that separates from the molecule, but instead of dissociating immediately it roams around the molecule for a considerable amount of time and extracts another hydrogen atom from the molecule prior to dissociation. This phenomenon has been coined roaming and has since been reported in the dissociation of a number of other molecules. In this paper we investigate roaming in Chesnavich’s \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {CH}_4^+$$\end{document}CH4+ model. During dissociation the free hydrogen must pass through three phase space bottleneck for the classical motion, that can be shown to exist due to unstable periodic orbits. None of these orbits is associated with saddle points of the potential energy surface and hence related to transition states in the usual sense. We explain how the intricate phase space geometry influences the shape and intersections of invariant manifolds that form separatrices, and establish the impact of these phase space structures on residence times and rotation numbers. Ultimately we use this knowledge to attribute the roaming phenomenon to particular heteroclinic intersections.


Roaming
For a long time it was believed that dissociation of molecules can only happen in two ways. Firstly, the original molecule can dissociate into smaller molecules and this is sometimes referred to as dissociation via the molecular channel. In order to dissociate, the system has to pass over a potential barrier representing the energy needed to break existing bonds and form new ones. Quantitative results on dissociation rates (or reaction rates in general) can be obtained via transition state theory. Alternatively an individual atom, called free radical, can escape from a molecule without forming new bonds and thus without passing over a potential barrier [1]. This is sometimes referred to as dissociation via the radical channel. Dissociation via both channels is well understood.
Recently, however, van Zee et al. [2] reported having experimentally observed dissociation of formaldehyde (H 2 CO) with CO in low rotational levels at an energy where dissociation through the molecular channel should have rather resulted in high rotational states of CO. The two proposed explanations for this behaviour are that either at least one of the vibrational modes of the transition state is quite anharmonic or there have to be two distinct molecular channels.
Townsend et al. [3] discovered in their study of formaldehyde (H 2 CO) a new form of dissociation that appears not to be associated with the molecular or the radical channel. In the process, an H atom separates from the molecule following the radical channel, but instead of dissociating it spends a considerable amount of time near the molecule and eventually abstracting the remaining H atom from the molecule to form H 2 . This type of dissociation is called roaming due to the nature of behaviour of the escaping H atom. No potential barrier or dynamical transition state is known to be involved in roaming.
The discovery of roaming stimulated extensive studies of formaldehyde photodissociation and roaming has since been accepted as the cause for the phenomenon observed by van Zee et al. [2].
Bowman and Shelper [4] have studied the dynamics of H 2 CO and CH 3 CHO to find evidence that roaming is more connected to the radical rather than the molecular channel. At the same time, roaming was observed at energies below the radical threshold.

Known results
In recent years dynamical systems theory has made an impact in chemistry by providing the means to understand the classical phase space structures that underlie reaction type dynamics [5][6][7]. This concerns the phase space geometry that governs transport across a saddle equilibrium point referred to as the molecular channel above. These ideas make precise the notion of a transition state which forms the basis of computing reaction rates from Transition State Theory which can be considered to be the most important approach to compute reaction rates in chemistry [8].
It is shown that surfaces of constant energy contain for energies above the saddle an invariant manifold with the topology of a sphere. This sphere is unstable. More precisely it is a normally hyperbolic invariant manifold (NHIM) which is of codimen-sion 2 in the energy surface. The NHIM can be identified with the transition state: it forms an unstable invariant subsystem located between reactants and products. What is crucial for the computation of reaction rates is that the NHIM is spanned by the hemispheres of another higher dimensional sphere which is of codimension 1 in the energy surface and which is referred to as a dividing surface. It divides the energy surface into a reactants region and a products region in such a way that trajectories extending from reactants to products have exactly one intersection with one of the hemispheres and trajectories extending from products to reactants have exactly one intersection with the other hemisphere. The construction of such a recrossing-free dividing surface is crucial for Transition State Theory where reaction rates are computed from the flux through a dividing surface which is computationally much cheaper than sampling trajectories.
For the global dynamics, the NHIM is significant because of its stable and unstable manifolds. The latter have the topology of spherical cylinders or 'tubes' which are of codimension 1 in the energy surface and hence have sufficient dimensionality to act as separatrices. In fact the stable and unstable manifolds separate the reactive trajectories from the non-reactive ones. The geometry of the stable and unstable manifolds and their location and intersections in the reactants and products regions carry the full information about the transition process including, e.g., state specific reactivity [9]. In the case of two degrees of freedom the NHIM is the Lyapunov periodic orbit associated with the saddle equilibrium point and the approach reduces the periodic orbit dividing surface (PODS) introduced earlier by Pechukas and Pollak and others [10,11].
Explaining the roaming phenomenon poses a new challenge to dynamical systems theory. The first attempts to use methods from dynamical systems theory related to ones underlying transition state theory to explain roaming can be found in the work of Mauguière et al. [12] in which they identified the region of the classical phase space where roaming occurs with the aid of numerous invariant phase space structures for Chesnavich's CH + 4 dissociation model [13]. They also introduced a classification of trajectories present in the system and matched them with the experimentally observed behaviour. The definition of roaming was formulated using the number of turning points in the radial direction in the roaming region. In light of a gap time analysis of Mauguière et al. [14], the definition was refined by means of the number of crossings of a dividing surface constructed in the roaming region. This refined definition is the best dynamical description of roaming to date.
Mauguière et al. [15] sudied the model of formaldehyde to find unstable periodic orbits in the roaming region. The homoclinic tangle of one such orbit was shown to be responsible for transport between two potential wells in a process that is closely linked to roaming. The periodic orbit involved do not arise as the Lyapunov periodic orbits associated with a saddle equilibrium points and the situation is hence different from the usual setting of transition state theory built on a potential saddle. Yet a recrossing free dividing surface can be constructed from such a periodic orbit. Such dividing surfaces may be other than spherical. It was shown in [16,17] that a spherical dividing surface near an index-1 critical point may bifurcate into a torus. In [18] a toric dividing surface was constructed near an index-2 critical point.
The authors of [19] found that the local geometry of the energy surface in an O 3 model may be toric and constructed a toric dividing surface using two unstable periodic orbits.
Even more recently, Huston et al. [20] report that they have found a correlation between the distribution of internal energies of CO and H 2 with the molecular channel and roaming. Particularly at lower energies, roaming trajectories have significantly more energy in H 2 . With increasing energy the differences decrease. Their definition of roaming is slightly different though, it involves rotation of H 2 around CO at a 'slowly varying and elongated distance'. The precise definition involves technical conditions on H 2 vibrational energy, time spent at a certain minimal distance from CO with low kinetic energy and large H-H bond length.

Objectives and outline of this paper
Because there is no single generally accepted definition of roaming, there is a clear need for a deeper understanding of the mechanisms behind dissociations.
In this work we present a detailed study of dissociation in the CH + 4 model by [13]. We discuss all types of dynamics present in this model and explain their connection to the underlying phase space geometry and invariant structures. We construct various surfaces of section and from the dynamics on these surfaces we deduce the role of invariant manifolds in slow dissociation and ultimately show a certain structure of heteroclinic tangles that causes roaming.
From the point of view of transition state theory we address two interesting problems. Firstly, it is not very well understood what happens in case reactants and products are divided by multiple transition states in series, which is a problem we address in this work. Secondly, we study the role of the local energy surface geometry in interactions of multiple transition states.
In the study we employ surfaces of section, all of which satisfy the Birkhoff condition [21] of being bounded by invariant manifolds. Using the surfaces of section we can observe dynamical behaviour such as roaming, but to understand the role of the local energy surface geometry and its implications to roaming, we generalise the Conley-McGehee representation [22][23][24] and study the dynamics on the energy surface in full 3 dimensions.
The paper is organized as follows. In Sect. 2 we introduce the Chesnavich's CH + 4 model. In Sect. 3 we discuss various periodic orbits and their role in setting up the problem of transport of phase space volumes between different phase space regions. In Sect. 4 we study the dynamics of the Chesnavich model by looking at trajectories from various perspectives. This section is followed by relating the dynamics to roaming in Sect. 5. The invariant manifolds that govern the dynamics and in particular roaming are discussed from a global perspective in Sect. 6. Conclusions are given in Sect. 7.
that is intended for the study of multiple transition states. In this model, only one H atom is free and the CH + 3 molecule is considered to be a rigid complex. It is a planar system and we study it in a centre of mass frame in polar coordinates (r, θ 1 , θ 2 , p r , p 1 , p 2 ), where r is the distance of the free H atom to the centre of mass (magnitude of the Jacobi vector), θ 1 is the angle between a fixed axis through the centre of mass and the Jacobi vector between CH + 3 and H, and θ 2 is the angle representing the orientation of CH + 3 with respect to the fixed axis. In these coordinates the kinetic energy has the form where m is the reduced mass of the system, and I is the moment of inertia of the rigid body CH + 3 . The system has a rotational SO(2) symmetry, which can be reduced giving a family of systems parametrised by the (conserved) angular momentum.
This reduction can be obtained from the following canonical transformation: Then p ψ = p 1 + p 2 =: λ is the total angular momentum and it is conserved. It follows that where U (r, θ) is the potential energy from [13] that we will discuss later in Sect. 2.3. In the last expression the term λ mr 2 p θ gives rise to a Coriolis force in the equations of motion.

General setting
As explained by [1], systems exhibiting roaming have a potential well for a small radius, representing the stable molecule, and with increasing distance between the dissociated components converges monotonously to a certain base energy, which we can assume to be 0. This is unlike the traditional bimolecular reactions that involve flux over a potential saddle. As shown in [19], under certain conditions the potential U admits two unstable periodic orbits that are not associated with any potential which, however, form the transition state to dissociation. We will find these orbits and use them to construct a toric dividing surface from them.
The argument for the existence of the periodic orbits is as follows. The dependence of the potential U (r, θ) on θ is due to the interaction between the anisotropic rigid molecule and the free atom. When r is sufficiently large, the potential U (r, θ) is essentially independent of θ . This is because for sufficiently large distances, the orientation of the CH + 3 molecule does essentially not influence the interaction with the free atom. Let us therefore assume for a moment that r is sufficiently large, so that U is rotationally symmetric and we can drop θ in the argument of U . The system reduced by the rotational symmetry then has the effective potential where p θ becomes a constant of motion. The reduced system admits a relative equilibrium, provided U (r ) is monotonous, U (r ) < 0 and U ∈ o(r −2 ) as r → ∞. Potentials of most chemical reactions, including CH + 4 → CH + 3 + H, meet this condition. The relative equilibrium is given by r = r p θ , p r = 0, where r p θ is the solution oḟ For the class of potentials U = −cr −(2+ ) with c, > 0, the relative equilibrium is unstable (Fig. 1). This follows from the reduced 1-degree-of-freedom Hamiltonian having a saddle at this equilibrium as can be seen from computing the Hessian matrix which is diagonal with the elements In the full system, the relative equilibrium is manifested as the unstable periodic orbits r = r p θ , p r = 0 and p ± θ such that Following general results on the persistence of normally hyperbolic invariant manifolds [25], these periodic orbit persist if the rotational symmetry is broken, provided the perturbation is not too big. Note that according to our assumptions these periodic orbits are not associated with a local maximum of U .
The condition U ∈ o(r −2 ) as r → ∞ is reminiscent of the assumption made by the authors of [26]. However, they consider a growth restriction near the origin, namely that for all θ λ 2 2mr 2 + U (r, θ) ∈ o(r −2 ) as r → 0, and additionally require λ 2 2mr 2 + U (r, θ) to have at most one maximum for each θ . We do not impose restrictions on U near r = 0 and admit several maxima.

Potential energy
The potential as suggested by Chesnavich [13] is the sum where U C H is a radial long range potential and U * a short range "hindered rotor" potential that represents the anisotropy of the rigid molecule CH + 3 [26,27]. The long range potential is defined by where x = r r e . The constants D e = 47 kcal/mol and r e = 1.1 Å represent the C-H dissociation energy and equilibrium bond length respectively. c 1 = 7.37 and c 2 = 1.61 result in a harmonic oscillator limit with stretching frequency 3000 cm −1 . A graph of U C H , using Chesnavich's choice of coefficients, can be found in Fig. 2. As expected for long range interactions, it is meant to dominate the potential for large values of r and not be subject to the orientation of CH + 3 . Therefore U C H is independent of the angle and its leading term for large r is r −4 . Since U C H also dominates the short range potential in the neighbourhood of r = 0, Chesnavich suggest a cut-off at r = 0.9. The cut-off is not near the region of interest in our study of roaming, nor does it have any significant implications.
The short range potential has the form where U 0 (r ) = U e e −a(r −r e ) 2 , is the rotor barrier, which is a smoothly decreasing function of the distance r , and U e = 55 kcal/mol is the barrier height, see Fig. 2. The constant a influences the value of r at which the transition from vibration to rotation occurs. The transition is referred to as early if it occurs at small r and as late otherwise. For comparison, see the late transition for a = 1, which we will be using, and the early transition for a = 5 in Fig. 3. Note the different proportions of the potential well (dark blue) with respect to the high potential islands along the vertical axis. Note that the angular dependence (1 − cos 2θ) in U * is π -periodic and even. These properties induce a reflection symmetry of U with respect to the x and y axes, because corresponds to the reflection about the x axis and corresponds to the reflection about the y axis.

Setting up the transport problem
We follow [12] and set a = 1 for a slow transition from vibration to rotation. In what follows we also assume λ = 0, unless stated otherwise. This section introduces features of the potential relevant to finding periodic orbits, defining dividing surfaces and formulating roaming in terms of transport between regions on the energy surface.

Energy levels and Hill regions
Here we give details about the features of the potential relevant to the dynamics of the system. Being the most basic characteristic of the potential, we look at critical points of the potential that give valuable information about local dynamics and at level sets that tell us about the accessible area in configuration space.
Due to the reflection symmetry of the potential about the x and y axes introduced above, critical points always come in pairs. We will denote them by q ± i , where i indicates the index of the critical point and the superscript + stands for the upper half plane θ ∈ [0, π), while − stands for the lower. Here we present a list of critical points: 27. The potential wells correspond to the two isomers of CH + 4 with the free H atom close to the CH + 3 molecule. All index-1 saddles are involved in isomerisation and the two index-2 saddles provide us with interesting geometries of the accessible regions in configuration space. For zero angular momentum (λ = 0), the critical phase space points of H are given by z ± i = (q ± i , 0) and z ± 1 = ( q ± 1 , 0). The critical energies are ordered as and all critical points can be found in the contour plot in Fig. 3.
For a given fixed energy E, we are interested in the accessible region in the configuration space which, following the celestial mechanics literature, we refer to as Hill region [28], and the geometry of the energy surface. Since the system as defined in Sect. 2.1 is a natural mechanical system H = T + U and the kinetic energy T is always non-negative, the Hill region consists of all (r, θ) such that U (r, θ) ≤ E and is bounded by the equipotential U (r, θ) = E.
To see what the Hill regions look like, we note that the two wells q ± 0 give rise to two topological discs that connect into an annulus at E = E 1 via q ± 1 . With E → 0 the annulus widens until at E = 0 it loses compactness and covers the whole plane except for a disc near the origin. This cut-out disc decomposes at E = E 1 into three, two areas of high potential around q ± Hill regions are shown for various energies in Fig. 3. For comparison, we also include Hill regions for the case a = 5 in Fig. 3, where the transition from vibration to rotation occurs earlier. Although energy levels remain topologically equivalent, note the larger potential well and the smaller energy interval where the boundary of Hill region consists of three circles.

Relevant periodic orbits
Next we study the invariant structures that can be found in the system at various energies. Critical points z ± i described in Sect. 3.1 are the most basic invariant structures at energies E i . In the following we discuss (non-degenerate) periodic orbits on the 3-dimensional energy surface. They create a backbone for the understanding the dynamical behaviour of our system.
Similarly to critical points, periodic orbits also come in pairs because of the symmetry of the potential. The periodic orbits are then related by the discrete rotational symmetry or the discrete reflection symmetry In contrast to critical points, non-degenerate periodic orbits persist in energy intervals forming one-parameter families. As periodic orbits evolve with varying energy, they occasionally bifurcate with other families of periodic orbits.
Based on the knowledge of Hill regions we gained in Sect. 3.1, we can formulate some expectations about periodic orbits in this system. For E ≤ E 1 , the system does not admit rotating periodic orbits, orbits that are periodic in θ and along which always p θ > 0 or p θ < 0. Rotating orbits project onto the configuration space as circles with the origin contained in their interior. Instead in the interval E 0 < E ≤ E 1 we can only expect vibrating, oscillator like, periodic orbits. Of special interest are periodic orbits that project onto a line with both ends on an equipotential. In the celestial mechanics literature these orbits are referred to as periodic brake orbits (the name is due to Ruiz [29]).
Here we present a list of the important families of periodic orbits together with a brief description of their evolution. Configuration space projections of the periodic orbits at E = 2 are shown in Fig. 4.
-Γ i : The family of periodic orbits Γ i is born in a saddle-centre bifurcation at energy E = −. 29. Until a host of bifurcations above E = 20, Γ i consists of hyperbolic brake orbits. Around E = 21.47 the orbits become inverse hyperbolic and at E = 22.27 they become heteroclinic to z ± 2 and undergo a Morse bifurcation, similar to those described in [16]. At higher energies, Γ i consists of rotating orbits that undergo further bifurcations. The periodic orbits are by some authors referred to as inner or tight periodic orbits. We denote the individual orbits by Γ i + and Γ i − .

Fig. 4 Configuration space projections of
, Γ a ± (red) and one orbit of the family Γ b (magenta) at energy E = 2 (Color figure online) ± is the brake orbit in the potential well associated with z ± 0 and for E > 22.27 the subscript ± corresponds to the sign of p θ along the rotating periodic orbit.
-Γ o : This family of unstable periodic orbits originates at r = ∞ at E = 0. With increasing energy the orbits monotonously decrease in radius and remain unstable until a bifurcation with Γ a and Γ b at E = 6.13, where Γ a and Γ b are described below. These periodic orbits are sometimes called outer or orbiting periodic orbits, because these are the periodic orbits with the largest radius at the energies where they exist. We denote the individual orbits with p θ > 0 and p θ < 0 by Γ o + and Γ o − respectively. -Γ a : These periodic orbits are created in a saddle-centre bifurcation at E = −.0602 as stable, turn unstable at E = −.009 and remain unstable until a period doubling bifurcation at 2.72. The family disappears in the aforementioned bifurcation with Γ o and Γ b . At all energies, the configuration space projection of Γ a is located between those of Γ i and Γ o and we will refer to the orbits as the middle periodic orbits. We denote the individual orbits with p θ > 0 and p θ < 0 by Γ a + and Γ a − respectively. -Γ b : The product of a saddle-centre bifurcation at E = −.0023 that quickly becomes inverse hyperbolic. Around E = 2.37 the orbits become elliptic and undergo a reverse period doubling at E = 2.4025. Note that the energetic gap between these two bifurcations is so small that they are almost indistinguishable in Fig. 5. After that Γ b remains stable until it collides with Γ a and Γ o . For E < 2.4025 the family consists of four periodic orbits with twice the period compared to all the previously mentioned ones. The orbits related by discrete symmetries mentioned above.
Γ i is important because its orbits lie in the potential well and have the largest radial coordinate r of all periodic orbits in the well. Γ o are the outermost periodic orbits and trajectories with a larger radial coordinate r and p r > 0 go to infinity in forward time, i.e. r → ∞ as t → ∞. As mentioned above, the configuration space projections of Γ a lie between Γ i and Γ o . In fact there are no other periodic orbits with single period (2π -periodic in θ ) in this region of configuration space. We use orbits of the family Γ a in Sect. 3.4 to define dividing surfaces and divide phase space into regions. Γ b is needed for a complete description of the evolution of Γ o and Γ a and its bifurcations may hint at qualitative changes of structures formed by invariant manifolds.
There are various other periodic orbits, most notably ones corresponding to stable vibrations of the bound CH + 4 molecule, Lyapunov orbits associated with z ± 1 and z ± 1 that play a role in isomerisation and periodic orbits involved in various bifurcations with the orbits mentioned above. All of these will not play a role in our further considerations.
With non-zero angular momentum, periodic orbits of a family remain related by the discrete rotational symmetry, but not by the discrete reflection symmetry and some other properties are different too. The inner periodic orbits are no longer brake orbits for λ = 0 and their projections onto configuration space are topological circles instead of lines. Similarly rotating orbits of the same family do not have the same configuration space projection and bifurcate at different energies. With increasing |λ| the differences become more pronounced.
In Fig. 5 we present the evolution of the orbits in above-mentioned families in the energy-action (E, S) plane. We will explain in Sect. 3.3 why flux through a dividing surface associated with a vibrating periodic orbit is equal to its action, while for rotating periodic orbits it is equal to twice its action. Figure 6 shows the evolution of the Greene residue of orbits in the families. The Greene residue, due to J. M. Greene [30] is a quantity characterizing the stability of the orbits. It is derived from the monodromy matrix, a matrix that describes the behaviour of solutions in the neighbourhood of a periodic orbit. The monodromy matrix and the Greene residue are defined as follows. For a periodic orbit Γ with the parametrisation γ (t) and period T , let M(t) be the matrix satisfying the variational equationṀ , with the initial condition The monodromy matrix is defined by M = M(T ). It describes how an initial deviation δ from γ (0) changes after a full period T . For δ sufficiently small the relationship is where Φ t H is the Hamiltonian flow. If δ is an initial displacement along the periodic orbit δ J ∇ H , then by periodicity δ is preserved after a full period T , i.e. Mδ = δ. A similar argument holds for an initial displacement perpendicular to the energy surface δ ∇ H . Consequently, two of the eigenvalues of M are λ 1 = λ 2 = 1. More details including a reduction of M can be found in [31].
As the variational equation satisfied by M(t) is Hamiltonian, the preservation of phase space volume following Liouville's theorem implies that the determinant det M(t) = det M(0) = 1 for all t. Therefore for the two remaining eigenvalues we have λ 3 λ 4 = 1 and we can write them as λ and 1 λ . Γ is hyperbolic if λ > 1, it is elliptic if |λ| = 1 and it is inverse hyperbolic if λ < −1.

Definition 1 The Greene residue of Γ is defined as
where M is the monodromy matrix corresponding to the periodic orbit Γ .

Transition states and dividing surfaces
In this section we discuss dividing surfaces associated with transition states, the backbone of Transition State Theory. Following [16,17] we define transition states more formally as follows.
Definition 2 (TS) A transition state for a Hamiltonian system is a closed, invariant, oriented, codimension-2 submanifold of the energy surface that can be spanned by two surfaces of unidirectional flux, whose union divides the energy surface into two components and has no local recrossings.
The name transition state is due to the fact it is a structure found between areas of qualitatively different types of motion, a transition between two types of motion so to say. One can imagine the transition between types of motion corresponding to physical states like reactants and products or the transition between rotation and vibration.
In a system with 2 degrees of freedom, a TS consists of unstable periodic orbit. Generally a TS is a codimension-2 normally hyperbolic invariant manifold, a manifold on the energy surface invariant under the Hamiltonian flow, such that instabilities transversal to it dominate the instabilities tangential to it [25,32].
In general, a dividing surface (DS) is a surface that divides the energy surface into two disjoint components. By a DS associated with a TS we mean a union of the two surfaces of unidirectional flux that is constructed as follows.
For a fixed energy E, let (r Γ , θ Γ ) be the projection of the periodic orbit Γ onto configuration space, then the DS is the surface (r Γ , θ Γ , p r , p θ ), where ( p r , p θ ) are given implicitly by the energy equation This construction also works for stable periodic orbits, but the resulting DS admits local recrossings. In the following a DS associated with a TS is always the surface constructed this way.
We will refer to the DSs associated with Γ i , Γ o and Γ a as inner, outer and middle, respectively. For our investigation, we do not need to distinguish between the DSs associated to Γ i + and Γ i − , therefore we always refer to the former unless explicitly stated otherwise. We are mainly interested in the influence of local energy surface geometry on the geometry of DSs and in the dynamics on DSs under the corresponding return map.
The geometry of the DSs is due to the form of the kinetic energy and the local geometry of the energy surface. It is well known that a DS associated to a brake periodic orbit is a sphere and the brake periodic orbit is an equator of this sphere [6]. The equator divides the sphere into hemispheres, whereby the flux through the two hemispheres is equal in size and opposite in direction. Trajectories passing this sphere from reactants to products intersect one hemisphere and the other hemisphere is crossed on the way from products to reactants. The flux through a hemisphere is then by Stokes' theorem equal to the action of the periodic orbit [23].
Rotating periodic orbits, on the other hand, such as Γ o , give rise to a DS that is a torus. The two orbits of the same family with opposite orientation are circles on the torus and divide it into two annuli with properties identical to the hemispheres. Using Stokes' theorem we find that the flux across each annulus is given by the sum of the actions of the two orbits, or simply twice the action of a single orbit [16].
Should it be necessary to distinguish the hemispheres or annuli of a DS by the direction of flux, the outward hemisphere or annulus is the one intersected by the prototypical dissociating trajectory defined by θ = 0, p r > 0, p θ = 0 and/or θ = π , p r > 0, p θ = 0. The inward hemisphere or annulus is then intersected by θ = 0, p r < 0, p θ = 0 and/or θ = π , p r < 0, p θ = 0

Division of energy surface
Using the inner and outer DSs we can define regions on the energy surface and formulate roaming as a transport problem.
The area bounded by the surface r = 0.9 and the two inner DSs represents the two isomers of CH + 4 . We denote the two regions by B + 1 and B − 1 . The unbounded region beyond the outer DS, denoted B 3 , represents the dissociated molecule.
It is therefore in the interaction region between the inner and the outer DS, denoted B 2 , where the transition between CH + 4 and CH + 3 +H occurs. When in B 2 , the H atom is no longer in the proximity of CH + 3 , but still bound to the CH + 3 core. This is the region, where the system exhibits roaming. Contained in B 2 are Γ a and various other periodic orbits that may play a role in roaming.
Dissociation can in this context be formulated as a problem of transport of energy surface volume from B 1 to B 3 . Such volume contains trajectories that originate in the potential well, pass through the interaction region and never return after crossing the outer DS. Since each trajectory passing from B 1 to B 2 crosses the inner DS and leaves B 2 by crossing the outer DS, we may restrict the problem to the interaction region. Because roaming is a particular form of dissociation, it too has to be subject to transport from the inner DS to the outer DS.
It is well known that transport to and from a neighbourhood of a unstable periodic orbits is governed by its stable and unstable invariant manifolds. The problem can be reformulated accordingly. This means, of course, by studying the structure of heteroclinic intersections of stable and unstable invariant manifolds of Γ i and Γ o , as well as with Γ a that, as we will soon see, sits inside the homoclinic tangle of Γ o .
We will denote the invariant manifolds of Γ i + by W Γ i + . We will further use a superscript s and u to label the stable and unstable invariant manifolds and add an extra superscript − and + for the branches that leave the neighbourhood of Γ i + to the CH + 4 side (r smaller) or to the CH + 3 +H side (r larger), respectively. W u+ Γ i + therefore denotes the unstable branch of the invariant manifolds of Γ i + that leaves the neighbourhood of Γ i + to the CH + 3 +H side. Invariant manifolds of other TSs will be denoted analogously. We remark that we may use TST to consider the evolution of periodic orbits in the energy-action plane shown in Fig. 5 in the context of transport of energy surface volume from B 1 to B 3 . Recall for Sect. 3.3 that the flux across the outer and middle DSs is twice the action of Γ o + and Γ a + respectively. The combined flux through both inner DSs is twice the action of Γ i + . We see that for E ≤ .32, the outer DS has the lowest flux, while for higher energies it is the inner DS.

Dynamics of the Chesnavich model
Before we proceed to the discussion of how invariant manifolds cause slow dissociations, let us describe some numerical observations of how the system behaves in certain phase space regions. The observations will later be explained using invariant manifolds. In the following, we offer insight into the amount of time needed to dissociate, the locations where dissociation is fast or slow and how these properties change with increasing energy. We use this knowledge to establish a link between invariant manifolds and slow dissociation on which we further elaborate in Sect. 5 in the context of roaming.

Residence times and rotation numbers
For various energies 0 < E < 6.13 where Γ o exists, we investigate trajectories starting in B + 1 , B 2 and B 3 on the surface θ = 0, p θ > 0. We study how long it takes trajectories to reach a terminal condition representing the dissociated state.
In Sect. 3.4 we said that we consider the molecule dissociated as soon as the system enters B 3 . Naturally, then the terminal condition should be that trajectories reach the outer DS. However using the outer DS raises uncertainty of whether a faster dissociation is a dynamical property or a result of the changing position of Γ o ± with energy. To prevent this uncertainty, we use a fixed terminal condition. Since for E → 0, the radius of Γ o ± diverges, no fixed terminal condition can represent the dissociated state for all energies. We decided to define the terminal condition by r = 15 that works well for E ≥ 0.4 at the cost of losing the energy interval E < 0.4.
In the following we consider residence times and rotation numbers, i.e. time and change in angle needed for trajectories starting on θ = 0, p θ > 0 to reach the surface r = 15 in B 3 .  Figure 7 shows rotation numbers for selected energies, with marked periodic orbits and invariant manifolds. As expected, initial conditions with p r > 0 large are the fastest ones to escape. The slowest ones are located near the periodic orbits and near p r = 0 (p θ large).
For E ≤ 2.5, almost all initial conditions in B + 1 were slow to escape. For higher energies, most of the slow dissociation occurs around Γ a + and Γ o + , the slowly dissociating trajectories have a negative initial p r very close to zero. This observation is easily explained by noting that configuration space projections of these trajectories are almost circular and spend most of the time in the region where the potential is very flat and almost independent of θ , thusṗ θ ≈ 0.
Chaotic structures that can be seen in B + 1 are the result of lengthy escape from a potential well. The only known structure responsible for fractal-like patterns and one closely linked to chaotic dynamics are invariant manifolds, in this case W Γ i + . Note that at E = 5, it seems that W Γ o + slows the dynamics down considerably more than W Γ i + . Rotation numbers, i.e. number of completed full rotations upon dissociation, closely match residence times suggesting that slowly dissociating trajectories are ones that rotate in B 2 and B 3 for a long time. More pronounced, due to the discrete nature of the number of rotations, are structures inside B + 1 , just below p r = 0 and in the neighbourhood of Γ o + and W Γ o + . Note in Fig. 8 that the fractal like structures recede with increasing energy and by E = 5 most of them lie either in B + 1 , near p r = 0 as mentioned above and in the proximity of the homoclinic tangle of Γ o + . The homoclinic tangle seems to tend to a homoclinic loop as it disappears for E → 6.13. It is also worth noting that fast and simple dissociation, i.e. low residence time and low rotation number, is not only becoming more dominant, but also speeding up, see Figs. 7 and 8 . Due to the increase in kinetic energy in the angular degree of freedom, the dissociating trajectories are naturally not becoming more direct with increasing energy.

Residence times on the inner DS
Similarly to the surface θ = 0, p θ > 0, we can study residence times and rotation numbers for trajectories starting on a DS. In Sect. 3.4 we formulated our problem as a transport problem from the inner to the outer DS. According to Sect. 3.3, trajectories enter B + 1 through one hemisphere of the inner DS and leave through the other. Naturally we are interested in the latter hemisphere.
Although it is not absolutely indispensable for qualitative purposes, we prefer to work on the DS in canonical coordinates. Due to the preservation of the canonical 2-form by the Hamiltonian flow, if we use canonical coordinates, the map from one surface of section to another is area preserving. Consequently areas of initial conditions on the inner DS corresponding slow or fast dissociation can be directly compared to the areas on the surface of section θ = 0, p θ > 0.
Canonical coordinates are obtained by defining a new radial variable ρ(r, θ) = r −r (θ ) that is constant along Γ i + , where the curver (θ ) is the approximation of the configuration space projection of Γ i + , similarly to [33]. Due to the symmetry of the system, Γ i + can be very well approximated by a quadratic polynomial for every energy. Next we use the generating function (type 2 in [34]) From that we obtain and therefore p σ = p θ +r (θ ) p ρ . The surface of section is now defined by ρ = 0, ρ > 0, i.e. the outward hemisphere of the inner DS corresponding to transport in the direction from B + 1 to B 2 . Figure 9 shows the distribution of residence times for initial conditions on the inner DS. We can see that slow dissociation is specific to two areas of the surface of section. Initial conditions on the rest of the surface leave B 2 quickly. Information from the two surfaces of section suggests that W Γ o + , and eventually W Γ a + , intersect the inner DS in the area with long dissociation times. The areas of slow dissociation are the most pronounced for low energies, at E = 2.5 they almost disappear. At E = 5 we see no sign of slow dissociation, the longest residence time found at the current resolution (6000×6000 initial conditions) was below 9. This suggests that the structure responsible for roaming disappears at an energy below 2.5. Note that even the slowest dissociation at E = 5 takes as long as the fastest ones at E = 1 or E = 2.
In summary we can say that the system exhibits various types of dissociation ranging from fast and direct, where the H atom escapes almost radially, to slow that involves H revolving a multitude of times around CH + 3 . Long dissociations seem to occur in fractal-like structures that are caused by invariant manifolds, proof of which will be given in Sect. 4.3.

Sections of manifolds
Let us now have a closer look at manifolds on the two surfaces of section presented above and establish a link between invariant structures and slow dissociation. In Sect. 4.1 we already noted that the homoclinic tangle of Γ i + is responsible for a fractal structure of slow dissociation of initial conditions in B + 1 . Furthermore the homoclinic tangle of Γ o + (and Γ o − ) is responsible for slow dissociation in the interaction region B 2 , especially at the top half of the energy interval. Residence times for initial conditions on the inner DS with outward direction for energies E = 1, 2, 2.5, 5. Note that the scale for E = 5 is different, because 9 is an upper bound for the residence time for initial conditions on the inner DS It is important to say that the section θ = const, p θ > 0 is not very well suited for the study of invariant manifolds. This is mainly due to the transition from vibration to rotation. The invariant manifolds W Γ i + may be nicely visible, but during this transition the invariant manifolds are not barriers to transport of surface area on this surface of section. Because parts W Γ i + rotate with p θ < 0 after the transition, they do not return to the surface of section. For the same reason there are trajectories that do not return to the surface of section. The return map associated with this surface of section is therefore not area preserving. This anomaly can be seen from odd shapes of invariant manifolds -heteroclinic points seem to be mapped to infinity.
Apart from the transition of W Γ i + from vibration to rotation, invariant manifolds may enter B ± 1 and be captured therein for a significant amount of time. Upon leaving B ± 1 the direction of rotation is unpredictable and this is true for invariant manifolds of all TSs. That is all we can say about the section θ = const, p θ > 0.
The section on the inner DS, just as all other DSs, does not suffer from these problems, because they do not depend on the direction of rotation. Moreover, these surfaces are almost everywhere transversal to the flow.   Fig. 10  and which is outside of it. For the majority of the area we note that σ = 0, p ρ > 0, p σ = 0 (equivalent to θ = 0, p r > 0, p θ = 0), the prototype of a fast dissociation, must lie inside W s− This problem is present on both the inner and outer DSs. Sections on both suffer from the fractal structure that is so characteristic for homoclinic and heteroclinic tangles. Heteroclinic points become visible after applying the return map at least once, the resulting tongues wind around the previously mentioned circles. On the downside, Γ a + is only hyperbolic until 2.72, therefore for higher energies the middle DS allows local recrossings. It can be still used as a surface of section and we can expect to see fewer heteroclinic points that cause tongues, but we need to keep local recrossings in mind.
In the next section we present a detailed view on the dynamics on the middle DS.

The observed dynamics and roaming
In this section we recall possible definitions of roaming used in previous works. We then elaborate on the observations above and analyse invariant manifolds on the middle DS with the aim to thoroughly explain how exactly roaming is linked to the heteroclinic tangles. Based on the explanation, a natural definition of roaming follows.

Roaming
Roaming in the chemistry literature refers to a kind of dissociation that is longer or more complicated than the usual dissociation with a monotonically increasing reaction coordinate that involves a saddle type equilibrium. While there is a sufficient amount of observations and intuitive understanding of what roaming is, an exact definition has not yet been generally adopted.
Mauguière et al. [12] proposed a classification of trajectories based on the number of turning points of trajectories in the interaction region B 2 . Later the authors refine their definition in [14] based on the number of intersections of a trajectory with the middle DS. Dissociating trajectories need to cross the middle DS at least three times before they are classified as roaming.
Huston et al. [20], on the other hand, set the criteria such that roaming trajectories have to spend a certain amount of time at a minimum radius, have low average kinetic energy and have on average a certain number of bonds over time.

The mechanism of roaming
Based on intersections of invariant manifolds, we would like to report on the types of trajectories in this dissociation problem and explain why the types exist. There is a general accord on the mechanism behind direct dissociation along the radical and molecular channel. The framework, that describes how codimension-1 invariant manifolds divide the energy surface in two and thereby separate reactive trajectories from non-reactive ones, is very well known in reaction dynamics, see [35][36][37][38].
Due to the different local geometries of the energy surface, we need to be careful with the invariant manifolds at this point. TSs that are brake orbits give rise to spherical DS and their invariant manifolds are spherical cylinders. TSs that are rotating orbits, just like ones belonging to the families Γ o and Γ a , give rise a toric DS that is based on two orbits instead of one. Therefore in the description of transport, invariant manifolds of both orbits have to make up a toric cylinder together. Invariant manifolds govern transport of energy surface volume as follows.
In CH + 4 → CH + 3 + H, we cannot discuss the molecular channel, but the radical channel and roaming is present. In general, if the H atom has enough kinetic energy to break bonds with CH + 3 , it escapes. Such a trajectory is contained in the interior of the invariant cylinder W u+ Γ i + , because it leaves the inner DS to the CH + 3 + H side. The same is true for W u+ Since the trajectory corresponding to θ = 0, p r > 0, p θ = 0 on the inner DS dissociates immediately, a part of W u+ All directly dissociating trajectories will be contained in the interior all of the above mentioned invariant cylinders. Moreover, directly dissociating trajectories are not contained in the interior of any other invariant cylinder.
As soon as a trajectory is contained in another cylinder, it is guided by that cylinder to cross the corresponding DS. Should a trajectory be contained in W s+ Γ i + , it will come back to the inner DS. In this way isomerisation, i.e. transport of energy surface volume between B + 1 and B − 1 , is possible via the intersection of the interiors of W u+ The intersections of the interiors of W s+ Γ a and W u+ Γ a or W s− Γ a and W u− Γ a , on the other hand, lead to the recrossing of the middle DS. In case a trajectory originating in B ± 1 dissociates after recrossing of the middle DS, by the definition of Mauguière et al. [14] it is a roaming trajectory. From the above it is clear that roaming trajectories are contained in intersection of the interiors of W u+ It remains to express the order of intersections of the DSs by a roaming trajectory with the invariant cylinders above.
In summary the arguments above enable us to say that, -directly dissociating trajectories are contained in W u+ and no other, -isomerisation and non-dissociating trajectories are contained in W u+ Note that since a trajectory contained in the cylinder W s− Γ a is automatically conveyed to W u+ Γ a after crossing the middle DS, we may omit mentioning one of the cylinders. A roaming trajectory could therefore be shortly characterized by W u+

Roaming on the middle DS
As mentioned in Sect. 4.3, the middle DS seems to be better suited for the study of roaming than the inner and outer DSs. More precisely, we will study dynamics on the outward annulus of the middle DS, i.e. the annulus crossed by the prototypical dissociating trajectory θ = 0, p r > 0, p θ = 0. We may introduce canonical coordinates on this annulus using a generating function in the same way as we did in Sect. 4.2, but for for the sake of simplicity we continue using the coordinates (θ, p θ ).
In the following elaboration we need means to precisely express the order in which invariant cylinders intersect the outward annulus of the middle DS. Based on the arguments in Sect. 5 The dynamics under the return map associated with the surface of section does not require W s+ Γ a for a complete and detailed description of dynamics. The simple fact that a point on the surface is mapped by the return map to another point on the surface is enough to deduce that the corresponding trajectory is contained W s+ Γ a and in fact, all other invariant cylinders made up of invariant manifolds of Γ a ± . Consequently, for a description of roaming on the outward annulus of the middle DS we only need W u+ By definition we have that contains trajectories that dissociate quickly. This is due to the fact that γ u+ i contains trajectories that just escaped from B + 1 and γ s− o contains those that reach the outer DS and therefore never return to the middle DS. Therefore points in γ u+ i ∩γ s− o do not have an image under the return map P, in fact the whole of γ s− o does not have an image. This is in accordance with the results on "reactive islands" by [35]. Figures 11, 12 and  13 show this intersection for various energies together with the first/last intersections of other invariant cylinders.
Note that trajectories passing through γ u+ i ∩ γ s− o reach the outer DS in varying amounts of time. We can expect the trajectory representing fast dissociation passing through θ = 0, p r > 0, p θ = 0 to take significantly less time than trajectories in the proximity of W s− Γ o , which may take arbitrarily long as they approach Γ o ± . Therefore if roaming was to be only defined by time spent in B 2 or in the neighbourhood of a periodic orbit, we can always find a suitable trajectory in γ u+ i ∩ γ s− o that is monotonous in r . Arguably, such a trajectory does not lead to an intramolecular hydrogen abstraction that has been reported in the context of roaming.
It remains to explain what happens to γ u+ Since trajectories corresponding to points in γ u+ i \ γ s− o do not dissociate, they return to the outward annulus of the middle DS unless they are asymptotic to a periodic orbit. The set γ u+ i \ γ s− o has an image under the return map P and it is Pγ u+ i . Note that the corresponding trajectories are guided to and from the middle DS by the invariant cylinders W s− Γ a , W u+ Γ a , W s+ Γ a , W u− Γ a . We remark that based on the understanding of lobe dynamics ( [38], trajectories that cross the outward annulus of the middle DS and do not reach the outer DS are repelled by Γ o ± and necessarily pass through the homoclinic tangle of Γ o ± . By the definition of Mauguière et al. [14], roaming trajectories cross the middle DS at least three times, which means crossing the outward annulus at least twice and the inward annulus at least once. Roaming trajectories must therefore contained in Pγ u+ i . In fact roaming trajectories that cross the outward annulus of the middle DS precisely n times before dissociating pass through Since recrossings of the middle DS are possible due to the homoclinic tangle of  In Fig. 13

Global study of the invariant manifolds that govern the dynamics
In this section we discuss an alternative way of studying dynamics on a 3-dimensional energy surface using the so called Conley-McGehee representation [22][23][24], described along with other alternatives in [39]. This is a very useful way of studying dynamics in full 3 dimensions, but to date has only been defined for subsets of energy surfaces that are locally a spherical shell. Since the Conley-McGehee representation does only works in B ± 1 of Chesnavich's CH + 4 model studied here, we introduce an extension of the Conley-McGehee representation that enables us to study energy surfaces with other geometry than in the Conley-McGehee case.

Conley-McGehee representation
The dynamics on the energy surface can be visualized in many ways. Just as it was done above, it can be viewed on various surfaces of section, most notably ones constructed around TSs. It is also possible to study the system locally using normal form approximations. The Williamson normal form [40] of the Hamiltonian in the neighbourhood of an index-1 critical point is for some λ, ω > 0. We found that for a fixed energy H 2 (q 1 , p 1 , q 2 , p 2 ) = h 2 , the energy surface can be locally viewed as a continuum of spheres parametrized by q 1 .
In the Conley-McGehee representation [22,24], is based on the spherical local geometry of an energy surface. While the normal form perspective above only applies locally, in the Conley-McGehee representation the whole energy surface is represented as a nested set of spheres parametrised in the radial direction by the reaction coordinate.
Advantages are immediate -the representation gives a global model of the energy surface and by construction reveals the spherical structure of the energy surface. For 2 degrees of freedom it enables us to study the 3-dimensional energy surface in the full 3 dimensions. Moreover, it enables to visualise the DSs as spheres that separate the energy surface into two disjoint components. It is also very natural that the flux through the hemispheres of the DSs is unidirectional and trajectories have to cross a particular hemisphere of the DS to pass from one component to the other.
Apart from DSs the Conley-McGehee representation enables us to visualise and therefore study TSs and their invariant manifolds that are spherical cylinders in a natural environment.

Toric extension of the Conley-McGehee representation
The Conley-McGehee representation in its original form applies to spherical geometries. The energy surface of the CH + 4 model, on the other hand, has a partially spherical and partially toric geometry, where many periodic orbits come in pairs and several DSs are tori. We therefore adapt the Conley-McGehee representation for the energy surface as follows.
The energy surface is defined by For very high energies, E > E 2 , where only r < 0.9 is energetically inaccessible due to the cut-off of the potential and E > U (r, θ) for all r ≥ 0.9, the whole energy surface has a toric local geometry. For any fixed radius r 0 and a fixed θ 0 , The radii of the concentric circles M E (r 0 , θ 0 ) on the ( p r , p θ )-plane depend on r 0 and θ 0 through U (r 0 , θ 0 ). The potential energy is not monotonous in r nor in θ . Recall from Sect. 3.1 that -the two wells q ± 0 are located at (1.1, 0) and (1.1, π) with U (q ± 0 ) = E 0 ≈ −47, -the two index-2 saddles q ± 2 are located at (1.63, π 2 ) and (1.63, 3π 2 ) with U (q ± 2 ) = E 2 ≈ 22.27.
Further recall from Sect. 2.2 that for r sufficiently large U (r, θ) is essentially independent of θ and r 2 U (r, θ) → 0 as r → ∞ for all θ .
It follows that the tori corresponding to r 0 = 1.1 and r 0 = r large , for some r large sufficiently large, intersect. This is because the circle M E (1.1, 0) has a larger radius than M E (r large , 0), while M E (1.1, π 2 ) has a smaller radius than M E (r large , π 2 ). The tori will always intersect if the radius is not a monotonous in r . In order to extend the Conley-McGehee representation, we need to reparametrise these tori so that their radii are monotonous in r for every θ . Define and P θ = U (r, θ)) p θ .

Now we have
The radius of the tori is monotonous in r and independent of θ and therefore the tori P 2 r + P 2 θ = r 2 foliating the energy surface are, unlike before the reparametrisation, disjoint in (θ, P r , P θ )-space.
Note that a section of the tori with a plane of section θ = θ 0 shows concentric circles, where the smallest one has the radius r = 0.9 due to the cut-off of the system. This is due to the fact that the boundary of the energy surface corresponds to a torus. Should it be desirable to have the whole (θ, P r , P θ )-space foliated by tori, it can be done by replacing r by r − 0.9 in the definitions of P r and P θ .

Extension to non-constant geometries
We remark that the construction above relies on the fact that at E > E 2 the energy surface has a purely toric geometry. We can slightly amend the construction to work for lower energies E ≤ E 2 , where the energy surface is not purely toric. E = U (r, θ) does not pose a problem for the definition of P r and P θ as it may seem on first sight. P r and P θ are only normalized conjugate momenta and by definition The momenta are therefore well defined on the whole energy surface.
Points on E = U (r, θ) are degenerate circles with radius 0 on the energy surface, but due to normalization correspond to circles in (P r , P θ ). Such a representation of the energy surface for lower energies is clearly flawed. For r large, we still have tori, but for smaller r , e.g. near Γ i , we do not see the spherical geometry we expect.
To solve this issue, we introduce different momenta in which the radius P 2 r + P 2 θ → 0 as U (r, θ) → E. We remark that the following only works for E ≤ E 1 . In the interval E 1 < E < E 2 , the projection of the energy surface on configuration space is the whole plane minus three discs, see Fig. 3. One is the potential energy cut-off and the other two are areas of high potential around index-2 critical points q ± 2 . Spherical and toric geometry cannot accurately represent a genus 3 surface.
Since for energies E < 0 the standard Conley-McGehee representation applies, we will restrict ourselves to the more interesting case 0 ≤ E ≤ E 1 . For the sake of simplicity, we retain the notation P r and P θ , making clear that we are discussing different momenta in a different energy interval than before. Let and P θ = r 3 1 2  U (r, θ)) for E = 0. The black circle corresponds to the radius r E defined in the text It follows that U (r, θ )).
While E −U (r, θ) makes sure that zero kinetic energy corresponds to P r = P θ = 0, the term r 6 seems perhaps less obvious. In the previous section we showed that it is important for P 2 r + P 2 θ to be monotonous in r for every θ . This is however not possible, because for some fixed angles θ = const the term (E − U (r, θ)) vanishes for several values of r . This is only possible in coordinates in which E = U (r, θ) correspond to coordinate lines for all E.
Let r E be the smallest r such that E = U (r, θ) has at most one solution for every θ for r E ≤ r . The term r 6 is the smallest even power of r such that r 6 (E − U (r, θ)) is monotonous in r on r E ≤ r for 0 ≤ E ≤ E 1 .
As a consequence of restricting the radius, the representation omits a significant part of B ± 1 . Since we formulated roaming as a transport problem from the inner DS to the outer DS, dynamics inside B ± 1 does not play a significant role in our study. All significant periodic orbits and DSs are well defined in the Conley-McGehee representation as presented here. Figure 14 shows the contour plot of r 6 (E − U (r, θ)) with a highlighted circle marking r 0 , the boundary of the representation defined above.

Consequences of the extension
In the representation as defined above, the geometry of the energy surface is preserved. For a fixed radius r 0 , we can see that the surfaces have the following topologies: -T 2 if U (r 0 , θ) < E for all θ , -a pinched torus if U (r 0 , θ) ≤ E for all θ and U (r 0 , θ 0 ) = E for some θ 0 , -S 2 ∪ S 2 if U (r 0 , θ 0 ) > E for some θ 0 . For E = 0, an example of each is shown in Fig. 15 in the canonical phase space coordinates (θ, p r , p θ ) and in the proposed extension of Conley-McGehee representation (θ, P r , P θ ). Note indeed in the latter that the surfaces are disjoint and present part of a foliation of the energy surface. We added an additional value of r = 4 to the extended Conley-McGehee representation to illustrate that the radius of the tori diverges for r → ∞, whereas it converges in the canonical phase space coordinates.
Due to the properties of the extended Conley-McGehee representation, we may study invariant structures and the aforementioned DSs globally on the energy surface. All techniques used to date either relied on surfaces of section or local approximations of the energy surface. In what follows, we study the structures on the energy surface in full three dimensions.
In Fig. 16 we present the periodic orbits Γ i + , Γ a ± and Γ o ± , all of which are TSs, and the associated DSs at energy E = 2.5. The construction of the DSs was discussed in Sect. 3.3. We remark that from a qualitative perspective Fig. 16 can be thought to represent the whole energy interval 0 < E < 2.72, where Γ a ± are unstable. The difference at higher energies is that Γ a ± is not a TS and the associated torus is not a DS.
We have used some of the DSs mentioned above in previous sections to study residence times, rotation numbers, and most importantly, the intersections of stable and unstable invariant manifolds of TSs. For the sake of clarity, in the following we left out the Γ i − and all associated structures, but everything said about Γ i + also holds for Γ i − . In the figures one can easily imagine another sphere just like the inner DS but shifted by π in the angular direction.
Note that here we take full advantage of the proposed extension of the Conley-McGehee representation to present structures on the energy surface with different local geometries -the inner DS is a sphere whereas the middle and the outer DSs are tori. This has to our knowledge not been done before. On the DSs we highlighted the respective TSs. Note that the inner DS is defined using one periodic orbit whereas the middle and the outer are defined using two. The individual orbits in the families Γ a and Γ o can be distinguished by the sign of P θ as they run in opposite directions.
The inner DS is divided by Γ i + into two hemispheres and the surfaces of unidirectional flux from the definition of a TS. Flux from B + 1 to B 2 crosses the hemisphere where predominantly P r > 0. The situation is similar for the annuli of the middle and outer DS. The two orbits divide the torus into two annuli, where the outward annulus is the one with larger P r .
The mechanism behind energy surface volume transport across DSs is governed by invariant manifolds of the corresponding TSs as discussed in the Sects. 4.3 and 5 . Here we present a new perspective for the study of invariant manifolds. In the following we use E = 5, because at this energy the TSs are evenly spaced and the lack of heteroclinic intersections facilitates understanding of this new perspective. Everything we say is applicable to the whole interval 0 < E < 6.13 relevant to roaming. Figure 17 displays the invariant manifolds of Γ i + and Γ o ± in full 3 dimensions. Note that the manifolds W Γ i + are only a sketch based on computed sections on the surface θ = 0. Computing the whole invariant manifold numerically is in this case relatively straight-forward, but for the sole purpose of illustration unnecessarily expensive.
Clearly visible is the structure of the manifolds, spherical cylinders formed by W Γ i One can clearly see that at E = 5, in fact in the whole energy interval E ≥ 2.5, the intersection of W u+ Γ i + with the middle and the outer DSs produces a topological circle centred at P θ = 0. This is purely the consequence of the spherical geometry induced by Γ i + . W Γ o ± for the same reason intersects the middle DS in two lines that should be seen as circles concentric with Γ a ± . This is in agreement with the sections Note that W u− Γ i + intersects the outward hemisphere of the middle DS at E = 2.5 in a shape that cannot be identified as a circle. This is due to the fact that we study W u− Γ i + that is asymptotic to Γ i in B + 1 as it leaves from B + 1 to B 2 . Moreover, it is deformed in the proximity of Γ a that exhibits a different kind of dynamics than Γ i , rotating as opposed to vibrating.
According to the findings in Sect. 5.3, the non-existence of roaming at higher energies due to the lack of intersection of W Γ i ± with W Γ o ± is immediate from Figs. 17 and 19.
The situation at E < 2.5 represented by Fig. 12 is very similar to the energy interval E ≥ 2.5. The main difference here are the heteroclinic intersections that cause roaming. These intersections are visible on the middle DS as well as in the extended Conley-McGehee representation. We remark that in the extended Conley-McGehee representation roaming occurs in the thin stripes between the two invariant cylinders W s− Γ o and W u− Γ o around P r = 0. Clearly the majority of the energy surface is occupied by more direct dynamics.

Conclusion
We have shown that numerical observations of long dissociation are caused by particular structures formed by invariant manifolds of TSs. These invariant manifolds are also responsible for multiple recrossings of the middle DS and consequently also for roaming.
We have shown that roaming trajectories that originate in the potential wells are captured in the homoclinic tangles of the outer TS and the middle TS dissociating. The transition for the potential wells into the homoclinic tangles is only possible in case the invariant manifolds of the inner and outer DS intersect and create a heteroclinic tangle. In case of Chesnavich's CH + 4 model, this heteroclinic intersection only exists for energies E ≤ 2.5 and therefore the system does not admit roaming at higher energy levels.
Our results can possible be directly extended to other chemical reactions, as the only significant assumption on the potential energy is that for all θ U (r, θ) ∈ o(r −2 ) as r → ∞.
This condition guarantees the existence of an outer TS thanks to which we may restrict roaming to a transport problem from potential wells representing a stable molecule to the DS associated with the outer TS.
Furthermore, our findings support then dynamical definition of roaming as suggested in [14].
We studied the above-mentioned invariant manifolds on various dividing surfaces, some of which highlighted the difficulties posed by the unusual energy surface geometry. This was the case even hough all surfaces of section satisfy the Birkhoff condition [21] of being bounded by invariant manifolds. These difficulties motivated us to construct an extension of the Conley-McGehee representation. Using the extension we were able to study the unusual energy surface geometry, TSs, the associated invariant manifolds and DSs in full three dimensions.