Phase space structures causing the reaction rate decrease in the collinear hydrogen exchange reaction

The collinear hydrogen exchange reaction is a paradigm system for understanding chemical reactions. It is the simplest imaginable atomic system with $2$ degrees of freedom modeling a chemical reaction, yet it exhibits behaviour that is still not well understood - the reaction rate decreases as a function of energy beyond a critical value. Using lobe dynamics we show how invariant manifolds of unstable periodic orbits guide trajectories in phase space. From the structure of the invariant manifolds we deduce that insufficient transfer of energy between the degrees of freedom causes a reaction rate decrease. In physical terms this corresponds to the free hydrogen atom repelling the whole molecule instead of only one atom from the molecule. We further derive upper and lower bounds of the reaction rate, which are desirable for practical reasons.


Introduction
We study the dynamics of the collinear hydrogen exchange reaction H 2 +H → H+H 2 , which is an invariant subsystem of the spatial hydrogen exchange reaction, using the potential provided by Porter and Karplus in [34]. In literature it is considered a paradigm system for understanding chemical reactions due to its simplicity and variety of exhibited dynamics. Because the system consists of three identical atoms confined to a line, it is the simplest imaginable system with 2 degrees of freedom modeling a chemical reaction.
The hydrogen atoms themselves are the simplest atoms in the universe. Because each consist of one proton and one electron only, an accurate potential energy surface for this reaction can be obtained via the Born-Oppenheimer approximation. Intriguingly enough, this system exhibits behaviour that is still not well understood.
The phenomenon we examine here is the counterintuitive observation that the reaction rate decreases as energy increases beyond a critical value. After all, one would expect to break bonds more easily using more energy. So far a satisfactory explanation of this phenomenon is missing and only an upper bound and a lower bound to the rate have been found. The upper bound is obtained by means of transition state theory (TST), due to [49]. TST is a standard tool for studying reaction rates due to its simplicity and accuracy for low energies, but it does not capture the decline of the reaction rate. The improvement brought by variational transition state theory (VTST) [13], does not capture this behaviour either.
Unified statistical theory, due to [23], which is in a certain sense an extension of TST to more complicated system, does capture the culmination of the reaction rate, but does not yield higher accuracy. The lower bound on the other hand does come quite close. It is obtained using the so-called simple-minded unified statistical theory [32].
A review of reaction rate results including TST can be found in [16]. Pechukas [28] and Truhlar and Garrett [40] review various extensions of TST.
Using lobe dynamics (introduced in [35]) we show how invariant manifolds of unstable periodic orbits guide trajectories in phase space. From the structure of the invariant manifolds we deduce that insufficient transfer of energy between the degrees of freedom causes a reaction rate decrease. In physical terms this corresponds to the free hydrogen atom repelling the whole molecule instead of only one atom from the molecule. We further derive bounds of the reaction rate, which are desirable for practical reasons.
In the remainder of this section we introduce the system, give an overview of TST and explain the current state of affairs with regards to the collinear hydrogen exchange reaction. Section 2 focuses on relevant periodic orbits and definition of regions of phase space. In Sect. 3 we introduce new coordinates using which we define a surface of section. In Sect. 4 we explain how we study invariant manifolds on the surface of section. In Sect. 5 we give a detailed insight into the structures formed by invariant manifolds and their role in the reaction. Section 6 is devoted to a novel way of breaking down heteroclinic tangles to provide a better understanding of the interplay of invariant manifolds of three TSs. In Sect. 7 we calculate various upper and lower bounds of the reaction rate.

Porter-Karplus potential
The collinear hydrogen exchange system consists of three hydrogen atoms confined to a line, as shown in Fig. 1, where r 1 and r 2 denote the distances in atomic units between neighbouring atoms. Forces between the atoms are given by the Porter and Karplus potential [34] is the standard potential for the hydrogen exchange reaction (collinear and spatial) used for example in [5,6,14,24,[31][32][33]. The system is considered to react, if it passes from the region of reactants (r 1 > r 2 ) to the region of products (r 1 < r 2 ) and remains there.
We point out two key properties of the Porter-Karplus potential: -the discrete reflection symmetry with respect to the line r 1 = r 2 , -the saddle point at r 1 = r 2 = R s := 1.70083.
The symmetry expresses the fact that we cannot distinguish between three identical hydrogen atoms, we can only measure distances between them. Hence, any statement referring to r 1 < r 2 automatically also holds for r 1 > r 2 .
Potential saddle points represent the activation energy needed for a reaction to be possible. In the all of this work we give energies as values in atomic units above the minimum of the system. In this convention the energy of the saddle point is 0.01456.
From a configuration space perspective, such a potential barrier is the sole structure separating reactants from products and the sole obstacle the system needs to overcome in order to react. This perspective implicitly assumes that the system does not recross the potential barrier back into reactants. Dynamical structures that cause recrossings are only visible from a phase space perspective. Figure 2 shows the potential energy surface near the potential saddle and cross sections of the potential at various values of r 2 . Due to diminishing forces between the atom and the molecule over large distances the differences between the cross sections fade after r 2 = 4 and are indistinguishable in double precision beyond r 2 = 40.

Definitions
The collinear hydrogen exchange reaction is described by the Hamiltonian H (r 1 , p r 1 , r 2 , p r 2 ) = p 2 r 1 + p 2 where p r 1 , p r 2 are the momenta conjugate to interatomic distances r 1 , r 2 , m H is the mass of a hydrogen atom and U is the Porter-Karplus potential described above. The equations of motion associated to H are as follows: The discrete symmetry of the potential translates into the invariance of H and the equations of motion under the map (r 1 , p r 1 , r 2 , p r 2 ) → (r 2 , p r 2 , r 1 , p r 1 ). The Hamiltonian flow generated by equations (2) preserves the energy of the system E = H (r 1 , p r 1 , r 2 , p r 2 ) and the phase space of this system is therefore foliated by energy surfaces H = E.
Examples of reactive and nonreactive trajectories are shown in Fig. 3. Note that nonreactive trajectories may cross the potential barrier in the sense that they cross the line r 1 = r 2 .
From the above it follows that the reaction rate at a fixed energy E can be calculated using a brute force Monte Carlo method as the proportion of initial conditions of reactive trajectories at infinity. Since the system decouples in a numerical sense around r 2 = 40, it is enough to sample a sufficiently remote surface in the reactants (r 1 > r 2 ) that is transversal to the flow, for example Since r 1 , r 2 is not a centre of mass frame, r 2 = const is not transversal to the flow. We remark that (r 2 , p r 2 − p r 1 2 ) are canonical coordinates on r 1 + r 2 2 = 50 that yield a uniform random distribution of initial conditions.

Transition state theory
Since its formulation in [49], TST became the standard tool for estimating rates of various processes not only in chemical reactions [16]. It has found use in many fields of physics and chemistry, such as celestial mechanics [10,15], plasma confinement [22] and fluid mechanics [25].
Key element of TST is the transition state (TS), a structure that is between reactants and products. There is no single generally accepted definition unfortunately, because in some publications concerning systems with 2 degrees of freedom TS refers to an unstable periodic orbit while in others TS is a dividing surface (DS) associated with the unstable periodic orbit. We adopt the following definition of a TS from [20]: 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 (the TS is the surfaces' boundary) of unidirectional flux, whose union divides the energy surface into two components and has no local recrossings.
For a system with 2 degrees of freedom as considered in this work, a closed, invariant, oriented, codimension-2 submanifold of the energy surface is a periodic orbit and it can be shown that the periodic orbit must be unstable [27,31,39]. In general, the TS has to be a normally hyperbolic invariant manifolds (NHIM), an invariant manifolds with linearised transversal instabilities that dominate the linearised tangential instabilities ( [8,11]).

Theorem 1 (TST) In a system that admits a TS and all trajectories that pass from reactants to products the DS precisely once, the flux across a DS is precisely the reaction rate.
We remark that in general the flux through a DS associated with a TS is an upper bound to the reaction rate [28,49].
Since its early applications, developments in the field led to a shift in the understanding of the TS to be an object in phase space rather than configuration space [41][42][43][44][45][46][47][48].
All relevant periodic orbits in this system are self-retracing orbits whose configuration space projections oscillate between equipotential lines, so called brake orbits ( [36]). As suggested by Pollak and Pechukas [31], let (r po 1 , r po 2 ) be the configuration space projection of a brake orbit at energy E, then the associated DS is the set of all phase space points (r For constructions of a DS near a saddle type equilibrium point in systems with more than 2 degrees of freedom see [41,47,48].
Hydrogen exchange results and evolution of understanding of TST follow.

Known results
In 1971, Morokuma and Karplus [24] evaluated three representatives of different classes of reactions. They found the collinear hydrogen exchange reaction to be the best suited for a study of the accuracy of TST due to smoothness, symmetry and simplicity. They found that TST agreed with Monte Carlo calculations up to a certain energy, but became inaccurate rather quickly after that.
In 1973 [29] Pechukas and McLafferty stated that for TST to be exact, every trajectory passing through the DS does so only once. In other words, TST fails in the presence of trajectories that oscillate between reactants and products.
In 1975 Chapman, Hornstein and Miller [5] present numerical results showing that transition state theory "fails substantially" for the hydrogen exchange reaction (collinear and spatial) above a certain threshold.
Pollak and Pechukas [31] proved in 1978 that flux through a DS constructed using an unstable brake orbit gives the best approximation of the reaction rate. In the presence of multiple TSs the authors introduce Variational TST (VTST)-using the DS with the lowest flux to approximate the reaction rate. These results detach TST from potential saddle points. The authors find for the collinear hydrogen exchange reaction that when TST breaks down, VTST can be significantly more accurate, even though both fail to capture the reaction rate decrease.
In 1979 Pollak and Pechukas [30] proved that TST is exact provided there is only one periodic orbit. Simultaneously, they derived the best estimate of the reaction rate so far for the collinear hydrogen exchange reaction in [32] using what they called Simple-minded unified statistical theory (SMUST).
Unified statistical theory (UST), due to Miller [23], attempts to take advantage of the difference of fluxes through all DSs and essentially treat regions of simple and complicated dynamics separately. The authors of [32] found that UST captures the drop in the reaction rate and elaborate on the deviation of UST from the actual rate. The 123 derivation of a lower bound (subject to assumptions) of the rate using the difference between TST and VTST is presented in the appendix of [32].
A rigorous lower bound is presented in [33]. It uses a DS constructed using a stable periodic orbit between two TS to estimate the error of TST. The accuracy of this lower bound for the hydrogen exchange reaction is remarkable.
In 1987 Davis [6] studied the hydrogen exchange reaction in phase space and considered the role of invariant structures. For low energies he showed that TST can be exact even if several TSs are present, provided that their invariant manifolds do not intersect. At higher energies he made some numerical observations of heteroclinic tangles of invariant manifolds and nearby dynamics. At high energies Davis found that a particular heteroclinic tangle grows in size and by assuming that it contains exclusively nonreactive trajectories he found a very accurate lower bound. The idea of this lower bound is very similar to [33], but Davis endures a computational cost to quantify trajectories instead of fluxes through DSs.
Davis also formulated an estimate of the reaction rate based on the observation that not many trajectories undergo a complicated evolution, as found by Pollak and Pechukas [32]. The estimate assumes that beyond a certain time dynamics in the heteroclinic tangle is randomised and 50% of the remaining trajectories are reactive.
Davis' observations hint at the crucial role played by invariant manifolds, but the precise manner in which this happens is not understood. Our aim is to explain the role of invariant manifolds in the reaction mechanism and extending it to the energy interval that Davis did not study, the interval with three TSs. We provide new understanding of the interactions between invariant manifolds of two and three TSs and consequently explain the counterintuitive reaction rate decrease.

Local geometry
Before we introduce periodic orbits that are relevant to the reaction mechanism, we describe the local energy surface geometry near a potential saddle point. We show that the neighbourhood necessarily contains an unstable periodic orbit and we highlight the importance of invariant manifolds to the local dynamics. The description remains true near unstable periodic orbits that do not lie near saddle points.
Consider the Williamson normal form [41,50] of a system near a saddle point. In the neighbourhood V of a potential saddle point, the system is accurately described in some suitable canonical coordinates (q 1 , p 1 , q 2 , p 2 ) by where λ, ω > 0. For a fixed energy H 2 = h 2 , this is equivalent to For a each fixed q 1 such that h 2 + 1 2 λq 2 1 > 0 this defines a sphere, as shown in Fig. 4. Depending on h 2 , the energy surface has the following characteristics: -If h 2 < 0, the energy surface consists of two regions locally disconnected near q 1 = 0, reactants (q 1 > 0) and products (q 1 < 0). -Reactants and products are connected by the saddle point for h 2 = 0. -For h 2 > 0, the energy surface is foliated by spheres. The radius of the spheres increases with |q 1 |. Locally the energy surface has a wide-narrow-wide geometry usually referred to as a bottleneck.
We remark that q 1 can be referred to as a reaction coordinate. To understand transport through a bottleneck, fix an energy h 2 slightly above 0 and consider the Hamiltonian equations for H 2 :q 1 = λ p 1 ,q 2 = ω p 2 , The degrees of freedom are decoupled with hyperbolic dynamics in (q 1 , p 1 ) and elliptic in (q 2 , p 2 ). Moreover q 1 = p 1 = 0 defines an unstable periodic orbit and q 1 = 0 defines a DS separating reactants from products. This DS, similarly to the one defined in Sect. 1.3, is a sphere that is due to the instability of q 1 = p 1 = 0 transversal to the flow and does not admit local recrossings. The sphere itself is divided by its equator q 1 = p 1 = 0 into two hemispheres with unidirectional flux -trajectories passing from reactants to products cross the hemisphere p 1 > 0, while trajectories from products to reactants cross p 1 < 0. Therefore q 1 = p 1 = 0 satisfies the definition of a TS. We remark that the DS can be perturbed and as long as its boundary remains fixed and transversality is not violated, the flux through the perturbed and unperturbed DS remains the same.
This description breaks down at high energies, when the periodic orbit may become stable, an event commonly referred to as loss of normal hyperbolicity. Then TST is inaccurate due to local recrossings of the DS. Loss of normal hyperbolicity occurs in the hydrogen exchange reaction, yet TST breaks down at lower energies due the presence of multiple transition states.
Having the same energy distribution between the degrees of freedom as the periodic orbit q 1 = p 1 = 0, its invariant manifolds are given by the stable being q 1 = −p 1 and the unstable q 1 = p 1 . They consist of two branches each-one on the reactant side with q 1 > 0, one on the product side with q 1 < 0. These manifolds are cylinders with the periodic orbit as its base. They are codimension-1 in the energy surface and separate reactive and nonreactive trajectories -reactive ones inside the cylinders and nonreactive outside Only reactive trajectories reach the DS. Note that in a configuration space projection, the separation between reactive and nonreactive trajectories is not as natural/obvious as in a phase space perspective. Therefore we study the structures made up of invariant manifolds that cause the reaction rate decrease in phase space.
We remark that bottlenecks are related to TSs rather than potential saddle points. Section 5 contains examples of bottlenecks unrelated to potential saddle points and a saddle point without a bottleneck.

Periodic orbits
For energies E above 0.01456, the energy of the saddle point, the system (1) admits periodic orbits that come in one-parameter families parametrised by energy. Initially we focus on each family separately and subsequently we investigate the interplay that governs the complicated dynamics exhibited by this system. We adopt the notation of [14] for different families of periodic orbits F n , where n ∈ N, and briefly describe their evolution with increasing energy. We remark that many families come in pairs related by symmetry and for simplicity we restrict ourselves to the r 1 ≥ r 2 half plane. We will refer to orbits of the family F n on the other half plane by F n . By F 0 we denote the family of Lyapunov orbits associated with the potential saddle, which as explained in Sect. 2.1 must be unstable for energies slightly above the saddle. The orbits lie on the axis of symmetry of the system r 1 = r 2 , see Fig. 5. Orbits of this family were used in TST calculations in many of the previous works.
A saddle-centre bifurcation at approximately 0.02204 results in the creation of two families-the unstable F 1 and the initially stable F 2 . The configuration space projections of these orbits are shown in Fig. 5. The unstable family F 1 is the furthest away from F 0 and does not undergo any further bifurcations. The F 2 family is initially stable, but undergoes a period doubling bifurcation at 0.02208 creating the double period families F 21 and F 22 . Unlike reported by Iñarrea and Palacián [14], we do not find these families disappear in an inverse period doubling bifurcation of F 2 at 0.02651. Instead F 21 and F 22 persist with double period until 0.02654, when they collide together with F 2 and F 0 , see Fig. 7. Consequently F 0 becomes stable. We would like to enhance the findings of [14] by remarking that F 21 and F 22 are briefly stable between switching from hyperbolic to inverse hyperbolic and vice versa, see Fig. 6.
At 0.02661, F 0 is involved in a bifurcation with a double period family F 4 that originates in a saddle-centre bifurcation at 0.02254. F 4 is a family symmetric with respect to r 1 = r 2 . For dynamical purposes we point out that above 0.02661 F 0 is inverse hyperbolic. Figures 6 and 7 show bifurcation diagrams of most of the families on the energyresidue and the energy-action plane. By residue R we mean the Greene residue as introduced by Greene in [9], where R < 0 means that the periodic orbit is hyperbolic, 0 < R < 1 means it is elliptic and R > 1 means it is inverse hyperbolic.
The residue is derived from a matrix that describes the local dynamics near a periodic orbit-the monodromy matrix. Let Γ be a periodic orbit with the parametrisation γ (t) and period T , and M(t) be the matrix satisfying the variational equation where The monodromy matrix is defined by M = M(T ) and it describes how a sufficiently small initial deviation δ from γ (0) changes after a full period T : where Φ t H is the Hamiltonian flow. According to Eckhardt and Wintgen, [7], if δ is an initial displacement along the periodic orbit δ J ∇ H , then δ is preserved after a full period T , i.e. Mδ = δ. Similarly an initial displacement perpendicular to the energy surface δ ∇H is preserved. Consequently, two of the eigenvalues of M are As (5) is Hamiltonian, the preservation of phase space volume following Liouville's theorem implies det M(t) = det M(0) = 1 for all t. Therefore the two remaining eigenvalues must satisfy λ 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.
where M is the monodromy matrix corresponding to the periodic orbit Γ .
Using (6) we can write R as Davis [6] mostly focused on the energy interval below 0.02214 and above 0.02655, the interval where TST is exact and the interval where two TSs exist, respectively.
In the light of normal form approximation described in Sect. 2.1, we remark that the approximation breaks down completely when F 0 loses normal hyperbolicity at 0.02655 at the latest. The loss of normal hyperbolicity is not the cause for the overestimation of the reaction rate by TST as it starts to deviate from the Monte Carlo rate well before 0.02300.

Phase space regions
We would like to give up the binary partitioning of an energy surface into reactants and products in favour of defining an interaction region inbetween into which trajectories can only enter once.
As explained in Sect. 2.1, TSs give rise to bottlenecks in phase space. Because F 1 gives rise to the bottleneck the furthest away from the potential barrier, we use it to delimit regions as follows. Denote DS 1 and DS 1 the DSs constructed using F 1 and F 1 according to Sect. 1.3. The interaction region is the region of the energy surface between the two DSs and it contains all other periodic orbits. Reactants and products are the regions on the r 1 > r 2 -side and the r 1 < r 2 -side of the interaction region respectively, see Fig. 8.
The advantages of this partition of space are immediate.
-All TSs and bottlenecks are in the interaction region or on its boundary. The dynamics in reactants and products has no influence on reactivity and to fully understand the hydrogen exchange reaction, it is enough to restrict the study to the interaction region. -Trajectories that leave the interaction region never return. This is true in forward and backward time. -It is impossible for a trajectory to enter reactants and products in the same time direction, unlike in the binary partitioning, where trajectories may oscillate between reactants and products.
123 Fig. 8 Regions in configuration space at energy 0.02400. The interaction region (red) bounded by two orbit from the family F 1 (blue), the region of reactants (blue) and the region of products (green). The orbit F 0 (black) is also included (Color figure online)

Definition of a Poincaré surface of section
Invariant manifolds are 2 dimensional objects on the 3 dimensional energy surface embedded in 4 dimensional phase space. To facilitate the study of intersections of invariant manifolds, we define a 2 dimensional surface of section on the energy surface that is transversal to the flow and intersects invariant manifolds in 1 dimensional curves.

Reaction coordinate and minimum energy path
Here we define a reaction coordinate, using which we can monitor the progress along a reaction pathway. Frequently a reaction coordinate is closely related to a minimum energy path (MEP) connecting the potential wells of reactants and products via the potential saddle. The coordinate as such is not a solution of the Hamiltonian system and, as remarked in [26], is of no dynamical significance to the system. A MEP can be defined as the union of two paths of steepest descend, the unique solutions of the gradient systeṁ one connecting the saddle (R s , R s ) to the potential well (∞, R min ), the other connecting (R s , R s ) to (R min , ∞). Figure 9 shows the MEP on a contour plot of U .

Surface of section
The MEP as defined above does not have an analytic expressing, but can be approximated using q 1 = 0, where as used by Davis [6] and shown in Fig. 9. Invariant manifolds are always transversal to the MEP and transversal to q 1 = 0 for the energy interval considered in this work. At higher energies Davis used q 1 = −0.04, q 1 = −0.07 and q 1 = −0.084 to avoid tangencies. We found that approximates the MEP significantly better, but a coordinate system involvingq 1 is rather challenging to work with. Throughout this work we use the surface of section Σ 0 defined by q 1 = 0,q 1 > 0. The conditionq 1 > 0 determines the sign of the momenta and guarantees that each point on Σ 0 corresponds to a unique trajectory. We remark that the boundary of Σ 0 does not consist of invariant manifolds and therefore it is not a surface of section in the sense of Birkhoff [4,Chapter 5].
For the sake of utility, we define the other coordinate q 2 such that (q 1 , q 2 ) is an orthogonal coordinate system on R 2 and the coordinate lines of q 2 are symmetric with respect to r 1 = r 2 . These conditions are satisfied by Note that q 2 = 0 is equivalent to r 1 = r 2 and q 2 is a reaction coordinate-it captures progress along q 1 = 0 and q 2 > 0 contains reactants, while q 2 < 0 contains products.

123
We remark that q 1 can locally considered a bath coordinate capturing oscillatory motion near the potential barrier. For a fixed energy, the energy surface is bounded in q 1 and unbounded in q 2 .

Symplectic coordinate transformation
Here we define a coordinate system in phase space, such that the coordinate transformation is symplectic. This requires finding the conjugate momenta p 1 , p 2 corresponding to q 1 , q 2 . For this purpose we use the following generating function (type 2 in [2]): Then One finds that From this we obtain This transformation has a singularity at r 1 = r 2 = R min , but U (R min , R min ) = 0.03845 is inaccessible at energies we consider. By straightforward calculation one finds that the symplectic 2-form ω 2 is indeed preserved: We remark that (q 2 , p 2 ) as defined above are the canonical coordinates on Σ 0 .

Transport and barriers
In this section we discuss the dynamics on the surface of section q 1 = 0 under the return map. This involves investigating structures formed by invariant manifolds via lobe dynamics due to [35].

Structures on the surface of section
The return map P associated with Σ 0 is defined as follows. Every point (q 0 , p 0 ) on Σ 0 is mapped to where T > 0 is the smallest for which q 1 (T ) = 0 along the solution with the initial condition where p 1 is given implicitly by the fixed energy E. P is symplectic because it preserves the canonical 2-form restricted to Σ 0 , see [3]. Because the Hamiltonian flow is reversible, P −1 is well defined. Each periodic orbit intersects Σ 0 in a single point that is a fixed point of P. Its stability follows from the eigenvalues of the monodromy matrix, as explained in Sect. 2.2. Due to conservation laws, the eigenvalues can be written as λ, 1 λ , 1, 1, see [7]. For TSs, the eigenvectors corresponding to λ, 1 λ define stable and unstable invariant manifolds under the linearisation of P near a fixed point.

Barriers formed by invariant manifolds
In the following we discuss invariant manifolds of TSs and their impact on dynamics with increasing energy. Let F i be a TS, we denote W F i its invariant manifolds as a whole, stable and unstable invariant manifolds are denoted W s F i and W u F i respectively. An additional +/− subscript indicates the branch of the invariant manifold with larger/smaller q 2 coordinate in the neighbourhood of F i , for example W s F i + and W s F i − . Recall from Sect. 2.1 that invariant manifolds of unstable brake orbits are cylinders of codimension-1 on the energy surface and they intersect Σ 0 in curves that divide Σ 0 into two disjoint parts each.
As mentioned in Sect. 2.2, the system has a single periodic orbit F 0 between 0.01456 and 0.02204. Its invariant manifolds do not intersect and act as separatrices or barriers between reactive and nonreactive trajectories, as shown at 0.01900 in Fig. 10. Reactive trajectories are characterised by a large | p 2 | momentum and are located above and below W F 0 . Nonreactive ones have a smaller | p 2 | momentum and are located between W s F 0 and W u F 0 . Consequently DS 0 , the DS associated with F 0 , has the no-return property and TST is exact ( [6]). F 1 and F 2 come into existence at 0.02204, but the reaction mechanism is governed entirely by W F 0 . W F 1 form a homoclinic tangle, but it only contains nonreactive trajectories. TST remains exact until 0.02215, when a heteroclinic intersection of W F 0 and W F 1 first appears. In the following we introduce the notation for homoclinic and heteroclinic tangles and subsequently introduce lobe dynamics due to [35] on the example of the homoclinic tangle formed by W F 1 , the F 1 tangle.

Definitions and notations
Let F i and F j be fixed points and assume W s F i and W u F j intersect transversally, as is the case in this system. The heteroclinic point Q ∈ W s F i ∩ W u F j converges to F i as t → ∞ and to F j as t → −∞. The images and preimages of Q under P are also heteroclinic points and therefore W s F i and W u F j intersect infinitely many times creating a heteroclinic tangle. If i = j, we speak of homoclinic points and homoclinic tangles.
Homoclinic and heteroclinic tangles are chaotic, since dynamics near its fixed points is locally conjugate to Smale's horseshoe dynamics (see [12] It should be clear that every tangle necessarily has pips. If Q is a pip, then Note that the end points of the segments are ordered, the first being closer to the fixed point along corresponding the manifold in terms of arclength on Σ 0 . Clearly P in pips other than the end points. Therefore P always maps lobes to lobes.

A partial barrier
Without knowing about invariant manifolds, the influence of a tangle on transport between regions of a Hamiltonian system may seem unpredictable and random. The role of invariant manifolds is well known and the transport mechanism may be intricate, yet understandable.
We explain this mechanism on the example of the F 1 tangle. The analogue in heteroclinic tangles will be apparent. The choice of the F 1 tangle at 0.02206 is due to the logical order in terms of increasing energy and its relative simplicity. Of the invariant manifolds, W s F 1 + and W u F 1 + form barriers similar to those discussed in Sect. 4.2 at all energies, while W s F 1 − and W u F 1 − form a homoclinic tangle. All branches of W F 1 lie in the region of nonreactive trajectories on the reactant side of F 0 , see Fig. 11.
Choose a pip Q 0 ∈ W s F 1 − ∩W u F 1 − , we will comment on the negligible consequences of choice later. The segments S[F 1 , Q 0 ] and U [F 1 , Q 0 ] delimit a region that we denote in reference to F 1 by R 1 . The complement to R 1 in the region bounded by W s F 0 + and W u F 0 + is denoted R 0 , see Fig. 12.
There is only one pip between Q 0 and P Q 0 , denote it Q 1 . In general the number of pips between Q 0 and P Q 0 is always odd (see [35]).
We define lobes using Q 0 , Q 1 and all of their (pre-)images. The way lobes guide trajectories in and out of regions can be seen on the lobe bounded by S[Q 1 , Q 0 ] and U [Q 0 , Q 1 ]. The lobe is located in R 0 , but its preimage bounded by S[P −1 Q 1 , P −1 Q 0 ] and U [P −1 Q 0 , P −1 Q 1 ] lies in R 1 . This area escapes from R 1 to R 0 after 0 iterations of the map P, we denote the lobe by L 1,0 (0). Analogously, by L 0,1 (0) we denote the lobe that is captured in R 1 from R 0 after 0 iterations and is bounded by S[P Q 0 , Q 1 ] and U [Q 1 , P Q 0 ]. We refer to images and preimages of L 1,0 (0) and L 0,1 (0) as escape lobes and capture lobes respectively. Note that due to the no-return property of the interaction region, escape and capture lobes cannot intersect beyond DS 1 .
Denote the lobe that leaves R i for R j , i = j, immediately after n iterations of the map P by L i, j (n).
In this notation we have for all k, n ∈ Z the relation Transition between R 0 and R 1 is closely connected to Q 0 and the transition from L i, j (1) to L i, j (0). All other lobes are confined by the barrier consisting of invariant manifolds to their respective regions. Near Q 0 , however, the barrier has a gap through which trajectories can pass. MacKay et al. [19] described this mechanism by saying that it "acts like a revolving door or turnstile." The term turnstile was born and lives on, see [21].
While W s F 1 − contracts exponentially near the F 1 , W u F 1 − stretches out. It is easy to see that S[F 1 , Q 0 ] is a rigid barrier-nearly linear and guiding all trajectories in its vicinity. W u F 1 − is a more flexible barrier in forward time -the manifold itself twists and stretches, alternately lying in R 0 and R 1 . The fluid shape of W u F 1 − is the result of complicated dynamics and the influence of S[F 1 , Q 0 ]. Stable manifolds behave similarly in backward time and the transition from rigid to flexible results in the turnstile mechanism.
The same is true for heteroclinic tangles. These imperfect barriers are responsible for nonreactive trajectories with high translational energy and reactive trajectories with surprisingly low translational energy. Due to this strangely selective mechanism we speak of a partial barrier.
Choosing any other pip than Q 0 for the definition of the regions merely affects the time in which lobes escape. Compared to definitions based on Q 0 , if we chose P Q 0 instead, escape/capture of lobes would be delayed by P, if we chose Q 1 , only escape lobes would be affected. This has implications for notation, not for dynamics or its understanding.

Properties of lobes
Here we state some of the basic properties of lobes that will be relevant in the following sections. The following statements assume that we study transport between two regions that are separated by a homoclinic tangle or a heteroclinic tangle and involves no other invariant manifolds. This provides useful insight into the complex dynamics of homoclinic and heteroclinic tangles.
If the intersection L i, j (0) ∩ L j,i (0) is non-empty, it does not leave the respective region and is not subject to transport. In this case we may redefine lobes to bẽ This justifies the following assumption.

Assumption 1
We assume that the lobes L i, j (0) and L j,i (0) are disjoint.
Equivalently we could assume L i, j (1) ⊂ R i and L i, j (0) ⊂ R j . In case of transport between several regions, we can only make statements based on the two regions that are separated by manifolds of the given tangle.
Each homoclinic and heteroclinic tangle involves a region bounded by segments of invariant manifolds, such as R 1 in Sect. 4.4. Since P is symplectic, almost all trajectories that enter the bounded region must eventually leave it. This can be formulated as Lemma 1 Let at least one of R i and R j be bounded. Then L i, j (0) can be partitioned, except for a set of measure zero O, as

Remark 1
The region R j has the no-return property iff escape lobes (L j,i ) are disjoint, or equivalently iff capture lobes are disjoint. Automatically then for all n > 0 Some of the intersections in Lemma 1 L i, j (0) ∩ L j,i (n) for n > 0 are empty sets. We are going to show that finitely many are empty at most.

Lemma 2 For all n
Using Fig. 12 as an example, Proof Without loss of generality assume R j is bounded and fix n 0 > 0. If We are going to argue that the only way for L i, j (−1) to reach L j,i (n 0 − 1) is by intersecting L j,i (n 0 ).
Denote Q 1 and Q 2 the pips that define L i, j (0) and P −n 0 Q 0 and P −n 0 Q 1 the pips L i, j (−1) lies inside R j (possibly partially in R i via another escape lobe) and so does U [P Q 1 , P Q 2 ], the part of ∂ L i, j (−1) that does not coincide with ∂ R j . Note that as all pips, P Q 1 , P Q 2 ∈ ∂ R j . The intersection point Q lies in the interior of the region bounded by and when mapped backward, Note for n 0 < 0, time reversal yields using a similar argument Following Lemmas 1 and 2, for k large enough L i, j (k) lies simultaneously in both regions forming a complicated structure. Since pips are mapped exclusively on ∂ R j , they aid identification of parts of lobes.
Due to (10), for n small we may study lobe intersections of the form that tend to be heavily distorted by the flow simply by mapping them forward or backward to less distorted intersections. However this does not work for for large k. On the other hand, we can expect the area of this intersection to shrink considerably with k, so their quantitative impact is limited. We remark that while almost the entire area of a capture lobe must escape at some point, this does not apply to entire regions. Regions may contain stable fixed points surrounded by KAM curves (sections of KAM tori) that never escape.
The picture of a heteroclinic tangle as a structure consisting of only two manifolds is oversimplified. In general heteroclinic tangles in a Hamiltonian system with 2 degrees of freedom can be expected to involve four branches of invariant manifolds. It takes four segments and two pips to define a region and consequently there will always be two turnstiles. The oversimplification is justified for tangles where the two turnstiles are made up of mutually disjoint lobes. Tangles with two intersecting turnstiles admit transport between non-neighbouring regions and we approach them differently.

Content of a lobe
In this section we use show how lobes guide trajectories in their interior.
Denote by μ the measure on Σ 0 , that is proportional to ω 2 Σ 0 (9). Under area preservation we understand that for any set A and for all k ∈ Z μ ( A) = μ P k A .
As a direct consequence of area preservation of a region we have for all k, n ∈ Z μ L i, j (n) = μ L j,i (k) . Assumption 2 Throughout this work we assume that μ(L i, j (0)) = 0.
Combining Assumptions 1 and 2 implies that L i, j (0), L j,i (1) ⊂ R j and if All other lobes may partially lie in both R i and R j , depending on the intersections of escape and capture lobes.

Definition 6
Assume R j is bounded. The shortest residence time in a tangle is a number k srt ∈ N, such that

Remark 2
The first lobe to lie partially outside R j is L i, j (−k srt ), because it intersects L j,i (0) ⊂ R i . The lobes L i, j (−k) and L j,i (k) are entirely contained in R j for 0 ≤ k < k srt .

123
Note that in a homoclinic tangle, since L i, j (−k) for 0 ≤ k < k srt must be mutually disjoint and all contained in R j , necessarily for all n > 0 and therefore L i, j (−k) intersects L j,i (0) ⊂ R i for all k > k srt . Due to reentries and Assumption 1, the statement is not true for L i, j (k) with k > 0, but an analogue holds in reverse time.
Reentries are possible in tangles where escape (and capture) lobes are not mutually disjoint, hence the following Lemma.
, P k 1 p ∈ R j and P k 3 −1 p ∈ R i follow from Assumption 1. Necessarily there exists k 2 , such that k 1 < k 2 < k 3 and p ∈ L j,i (k 2 ). Since k 2 may be different for every p, the union over k 2 follows.
The argument can be easily generalised for tangles that govern transport between multiple regions. One only needs to observe that p can return to R i from any region.
In the F 1 tangle at 0.02215, reentries can be deduced from the intersection L 0,1 (1)∩ L 1,0 (0) that lies completely in R 0 . See Fig. 13 for comparison of a tangle at 0.02215 with reentries and at 0.02210 without. Note that both tangles have k srt = 1.
Instantaneous transport between regions is described by the turnstile mechanism. Transport on a larger time scale can be studied using a measureless and weightless entity (species, passive scalars or contaminants [37,38]) that is initially contained and uniformly distributed in a region, as done in [35]. Its role is to retain information about the initial state without influencing dynamics indicate escapes and reentries via lobes.
The challenge of studying lobes over large timescales is to determine which regions a lobe lies in and correctly identifying the interior of a lobe. For this we propose a partitioning of heteroclinic tangles into regions of no return outside of which the evolution of lobes is of no interest.

Influence of tangles on the reaction rate
In this section we discuss the evolution of homoclinic and heteroclinic tangles in the entire energy interval 0 < E ≤ 0.03000 and their influence on dynamics in the interaction region. The dynamics for higher energies is due to the lack of bifurcations analogous. The study of invariant manifolds employs lobe dynamics and a new partitioning based on dynamical properties. An in-depth review of invariant manifolds in Fig. 13 The F 1 tangle at 0.02210 (above) and at 0.02215 (below). Both homoclinic tangles have k srt = 1, that can be seen by L 0,1 (0) ∩ L 1,0 (1) = ∅ shown in cyan. At 0.02215 the tangle admits reentries a chemical system and structural changes in tangles caused by bifurcations has to our knowledge not been done before.

Energy interval where TST is exact
TST is exact in the presence of a single TS (due to [30]) and remains exact in case of multiple TSs provided their invariant manifolds do not intersect (due to [6]). Therefore results of TST and Monte Carlo agree on the interval from 0 to 0.02215. W F 0 separate reactive and nonreactive trajectories, see Sect. 4.2, while the F 1 tangle captures nonreactive trajectories only. Some properties of the F 1 tangle are carried over to higher energies, such as shape of lobes or k srt . Figure 14 shows W F 0 and W F 1 approaching prior to the intersection at 0.02215 and the failure of TST.
Each change of structure seems to coincide with a bifurcation of a periodic orbit. The decrease k srt from 3 to 1 over the energy interval, shown in Fig. 15, coincides with the period doubling of F 2 at 0.02208 and the period doubling of F 21 before 0.02209. From a quantitative perspective, the tangle and its lobes grow larger in area.

Point where TST fails
At 0.02215, W F 0 and W F 1 interact through heteroclinic intersections. Instead of minor changes in the overall topology of the invariant manifolds, we come across something that is better described as a chain reaction.
Firstly, we observe that W F 0 + and W F 1 − intersect forming a heteroclinic tangle, see Fig. 16. Consequently, TST starts to fail (see [6]) and the Monte Carlo reaction rate is lower than TST. W s F 0 and W u F 0 form a partial barrier and this enables the F 1 tangle to capture reactive trajectories. We also find heteroclinic intersections of W F 1 − and W F 1 + as shown in Fig. 17, as well as W F 1 − and W F 0 − . Recall that statements for F 1 also hold for F 1 .
Choose two pips in the F 0 -F 1 tangle, so that the region bounded by W F 0 + and W F 1 − denoted R 0 satisfies R 1 ⊂ R 0 (Fig. 16) and define R 0 using symmetry. As L 0,1 (0) and L 1,0 (1) in the F 1 tangle contain heteroclinic points that converge towards F 0 (forward or backward time), they necessarily intersect in R 0 (see Fig. 13). By definition, L 0,1 (0) ∩ L 1,0 (1) contains trajectories that reenter R 1 after they have escaped and consequently R 1 (and R 1 ) loses its no-return property. In particular, trajectories that periodically reenter R 1 may exist and if they do, they will be located in L 0,1 (0) ∩ L 1,0 (k) ∩ . . . for somek.
By symmetry L0 ,1 (0) and L0 ,1 (1) also contain heteroclinic points that converge towards F 0 and they cannot avoid intersecting L 1,0 (1) and L 0,1 (0) respectively. Figure 17 portraits the intersecting invariant manifolds. These intersections guide trajectories that may cross DS 0 multiple times and result in an overestimation of the reaction rate by TST. Due to the size of the lobe intersections, the overestimation is small but increases with energy. VTST suffers from recrossings too as it estimates the rate using the DS with lowest flux, but none of the DSs is recrossing-free.
Due to a high k srt and small area of lobes, we avoid details of the F 1 -F 1 tangle until higher energies. We remark that lobes in the F 1 -F 1 tangle do not intersect outside of the bounded region.

Definitions of important regions
We have established that TST fails at 0.02215 due to recrossings. In this section we give a detailed description of homoclinic and heteroclinic tangles at 0.02230 and explain the transport mechanism in these tangles using lobes. The energy 0.02230 is representative for the interval between TST failure at 0.02215 and one of several period doubling bifurcations of F 21 at 0.02232. Moreover, lobes at 0.02230 are sufficiently large to study.
For the sake of simple notation, in what follows Q 0 , Q 1 , Q 2 and Q 3 denote pips that differ from tangle to tangle. To avoid confusion, we always clearly state which tangle is discussed.
First we discuss the homoclinic tangles of F 0 , F 1 and F 1 at 0.02230. We define regions relevant to these homoclinic tangles shown in Fig. 18 as follows.
Denote R 0 , the region bounded by W F 0 + and W F 1 − . The F 0 -F 1 tangle is responsible for most of the complicated evolution of reactive trajectories at 0.02230. The regions above and below the F 0 -F 1 tangle are R 2 and R 3 respectively.
The region inside the F 1 tangle bounded by W F 1 − is denoted R 1 . Further we denote R 4 the region bounded by W F 0 + that is relevant for the F 0 tangle. A near-intersection of W F 0 + in R 1 suggests that R 4 is smaller after the period doubling bifurcation of F 21 at 0.02232.

Homoclinic tangles
First we concentrate on the F 0 tangle at 0.02230, followed by the F 1 tangle, both depicted in Fig. 19. In both it is possible to identify a number of lobes that explain the dynamics within.
By far the largest intersection in the F 0 tangle is L 3,4 (−1) ∩ L 4,2 (2). It comprises most of the white area in R 4 occupied by nonreactive trajectories and we can deduce the structure of the intersection from L 3,4 (0) and L 4,2 (1) as follows. As an image of L 3,4 (0), the larger part of L 3,4 (−1) is bounded by S[P Q 1 , P Q 0 ] ∪ U [P Q 0 , P Q 1 ] with pips indicated in Fig. 19. This is nearly a third of the entire region Thanks to pips we are able to deduce that the majority of trajectories in the F 0 tangle is due to the intersection of these two lobes.
Note that part of an escape lobe extends to the product side of F 0 and contains reactive trajectories. This part of the lobe enters R 4 via L 3,4 (1), most of which is mapped to L 3,4 (0) ∩ L 4,2 (2) and escapes into R 2 via L 4,2 (1). Using an analogous argument we find that the part of a capture lobe lies on the product side of F 0 and carries reactive trajectories that escaped from R 4 .
The F 1 tangle has only one pip between Q 0 and P Q 0 and therefore a simpler structure. L 0,1 (0) ∩ L 1,0 (1) implies k srt = 1, therefore trajectories pass through this tangle quickly. Most nonreactive trajectories of the F 0 tangle pass inbetween L 1,0 (0) and L 0,1 (1) and avoid the F 1 tangle. This follows from its adjacency to Q 0 , which is only mapped along the boundary of R 1 always on the reactant side of F 0 . Similarly we can follow the area between L 1,0 (0) and L 0,1 (2) on the product side of F 0 using the F 1 tangle and symmetry.
The considerable size of lobes on the product side of F 0 carries information about nonreactive trajectories. The part of L 0,1 (1) on the product side of F 0 enters R 1 via the upper part of L 0,1 (0), just above the indicated intersection with L 1,0 (−1). Since this area does not lie in L 1,0 (1), it is has to be mapped to L 0,1 (−1)\L 1,0 (0) that remains in R 1 and is defined by the pips P Q 1 and P 2 Q 0 located on S[F 1 , P Q 0 ]. Further this area will be mapped in L 1,0 (1)\L 0,1 (0) and, unlike the part of L 1,0 (1) bordering S[P −1 Q 1 , P −1 Q 0 ], back into products.
As energy increases, we observe that the nonreactive mechanism of the F 0 tangle grows slower than the nonreactive mechanism in the F 1 tangle or even shrinks. The Fig. 19 Homoclinic tangles associated with F 0 and F 1 respectively at 0.02230 later involves crossing the axis q 2 = 0, which on Σ 0 coincides DS 0 . Due to symmetry the same happens in the F 1 tangle. Therefore the flux across DS 0 grows twice as quickly as across DS 1 . Therefore eventually DS 1 becomes the surface of minimal flux.

Heteroclinic tangles
Heteroclinic tangles partially share shapes, lobes and boundaries with homoclinic tangles and their description of transport must agree. Recall heteroclinic tangles have two turnstiles and two sets of escape and capture lobes.
For the sake of simplicity, we rely on pips and prior knowledge from Sect. 5.4 to interpret Fig. 20. Define R 0 in the F 0 -F 1 tangle using W F 0 + and W F 1 − and the pips 123 Fig. 20 The F 0 -F 1 tangle and the outline of F 1 -F 1 tangle at 0.02230 Q 0 and Q 2 . A single pip is located on ∂ R 0 between Q 0 and its image, the same is true for Q 2 .
It is worth mentioning that the lobes governing transport from R 2 to R 3 , L 0,3 (1) and L 2,0 (0), are disjoint. Nonreactive trajectories originating in R 2 spend some time in R 0 . This agrees with our conclusions on the nonreactive mechanism in the F 1 tangle.
The reactive mechanism in the F 0 -F 1 tangle involves the capture lobe L 3,0 (1) part of which is mapped to L 3,0 (0) \ L 0,2 (1) and on to L 3,0 (1) ∩ R 0 , part of which lies in L 0,3 (1). The area of this intersection is small in R 0 .
Understanding the F 1 -F 1 tangle is very involved, as the boundary of the tangle requires several segments of W F 1 − and W F 1 + . We propose a different point of view. In all tangles above, we have found that escape from the bounded region in a tangle, all area above the uppermost and below the lowermost stable invariant manifold escapes without further delay. For example in the F 0 -F 1 tangle, L 0,2 (1) located above W s F 1 − and L 0,3 (1) located below W s F 0 + escape to reactants and products respectively, because as the stable manifold bounding the lobe contracts, the unstable manifold is unobstructed to leave the interaction region. In this sense that we propose only stable invariant manifolds to be considered a barrier in forward time.
Using this reasoning, concentrate on the area between S[ F 1 , Q 0 ] and S[F 1 , Q 0 ] in the F 1 -F 1 tangle. Everything above S[Q 3 , Q 2 ] and below S[ Q 3 , Q 2 ] may pass through the tangle, but evolves in a regular and predictable manner from R 3 to R 2 or vice versa. We remark that this area is the intersection of two turnstiles. The same argument applies to the areas above S[Q 1 , Q 0 ] and below S[ Q 1 , Q 0 ]. Complicated dynamics is restricted to R 1 , as defined in the F 1 tangle, R 1 and an island near F 0 and should be treated separately from predictable areas.
Using this line of thought enables us to formulate bounds and estimates of the reaction rate. Before we proceed to quantitative results, we conclude this section by describing the evolution of tangles with increasing energy.

Higher energies
Based on the analysis in Sects. 5.4 and 5.5 for tangles at 0.02230, here we discuss the evolution of tangles at higher energies and their impact on dynamics in the interaction region. As the mechanisms have been described, most of our comments concern sizes of lobes and duration of escape from a tangle.
An interesting question arises from the connection between bifurcations and changes in geometry of invariant structures. The causal relationship is not evident. Also bifurcations are mostly thought of as local events. However as they seem to affect invariant manifolds, a change in tangles propagates instantaneously throughout the whole space. This phenomenon reminds of the infinite propagation speed in the heat equation.
The next bifurcation above 0.02230 according to Sect. 2.2 is a period doubling of F 21 at 0.02232, followed by a saddle-centre bifurcation that creates F 3 and F 4 at 0.02254 and a bifurcation of F 3 where F 31 and F 32 are created at 0.02257. At around 0.02523 follows another period doubling of F 21 , F 21 collides with F 2 at 0.02651 and subsequently F 2 collides with F 0 at 0.02654.
The major consequence of the bifurcation of F 21 at 0.02232 is a new intersection of W u F 0 + and W s F 0 + labeled Q 0 in Fig. 21. This reduces the number of pips between Q 0 and P Q 0 to one and therefore lobes are no longer made up of two disjoint sets. The F 0 tangle resembles the F 0 -F 1 tangle at 0.02230. Also the size of R 4 is reduced.
123 Fig. 21 The F 0 tangle and the F 1 tangle at 0.02253 In the F 1 tangle we see L 0,1 (2) cross DS 0 twice as shown in Fig. 21. All L 0,1 (k) for k > 2 and also L 1,0 (k) with k < −2 therefore pass through R 1 . Moreover, the tip of L 0,1 (2) approaching R 1 can be expected to pass cross R 1 after the bifurcations at 0.02254 and 0.02257. A small remark regarding notation. At this energy L 0,1 (2) lies in R 0 , R 0 , R 1 , R 1 , R 2 and R 3 , but we maintain the notation for consistency.
At 0.02400, L 0,1 (2) in the F 1 tangle passes through R 1 twice and the number increases at higher energies. Almost all lobes lie in almost all regions, but the mechanism for fast entry and exit of the tangles remain the same. Figure 22 shows R 0 and R 0 . While R 4 is considerably larger than R 1 at 0.02253, the opposite is true at 0.02400. Recall that R 4 contains predominantly nonreactive trajectories that do not cross DS 0 , Fig. 22 Indication of boundaries of R 0 and R 4 (above) and the F 1 tangle at 0.02400 (below). The area in the F 1 tangle highlighted in cyan is the part of L 0,1 (0) ∩ L 1,0 (1) that originates in products and is guided by W s whereas R 1 mostly contains ones that do. The overestimation of the reaction rate follows.
The capture lobes in the F 1 tangle guide predominantly trajectories from products into R 1 , as shown in Fig. 22. A significant portion of R 1 is taken up by L 0,1 (0)∩L 1,0 (1) and it is prevented by W s F 1 − from escaping into reactants. Moreover, the a large part of the intersection lies below W s F 1 + , see Fig. 22, that guides it into back products as W s F 1 + contracts. Heteroclinic tangles mirror the changes of the homoclinc tangles (Fig. 23).

Loss of normal hyperbolicity
F 0 loses normal hyperbolicity and becomes stable at 0.02654, in a bifurcation involving F 2 , F 2 , F 21 , F 21 . TST cannot be based on F 0 and W F 0 cease to exist. The sudden disappearance of invariant manifolds has no dramatic consequences. As can be deduced from Fig. 23, W F 0 are at energies below 0.02654, very close to W F 1 − and W F 1 + and naturally take over the role of W F 0 . Throughout the energy interval from 0.02206 when F 1 appears to the loss of normal hyperbolicity at 0.02654, we see a transition of dominance from F 0 to F 1 -F 1 .
The loss of normal hyperbolicity of F 0 simplifies dynamics due to the presence of fewer TSs, for example compare Figs. 23 and 24.
At 0.02661, F 0 collides with F 4 and becomes inverse hyperbolic. Due to the inverse hyperbolicity, W F 0 exist, but they must contain a twist that is manifested as a reflection across the F 0 (see [1]), i.e. have the geometry of a Möbius strip. At the same time W F 0 are enclosed between W F 1 − along with W F 1 + , but with cylindrical structure. Consequences of the geometry of W F 0 are unknown.
There are no more significant bifurcations above 0.02661 and therefore apart from growing tangles and lobes, the tangles remains structurally the same.
Together with W F 0 we observe the disappearance of R 4 and of the mechanism that carries nonreactive trajectories through the F 0 tangle without crossing DS 0 . Consequently, all trajectories that pass through the F 1 − F 1 tangle cross DS 0 at least twice. Each hemisphere of DS 1 still possesses the no-return property, which means trajectories cross DS 1 at most twice. Trajectories that avoid the tangle cross both DSs once or not at all.
Similarly to lower energies, R 1 is predominantly made up of L 0,1 (0) ∩ L 1,0 (1) in the F 1 tangle or L 2,0 (0) ∩ L 0,3 (1) in the F 1 -F 1 tangle, as shown in Fig. 24. The argument that trajectories in the F 1 -F 1 tangle below and above all stable manifolds leave the interaction region is still valid. Capture lobes are disjoint, therefore it is not possible to reenter the bounded region. Although R 1 and R 1 admit return, R 1 ∪ R 1 possesses the no-return property.

Known estimate
Davis [6] formulated bounds and an estimate of the reaction rate based on numerical observation of dynamics. He observed that a significant portion of trajectories leave the heteroclinic tangle above 0.02654 after one iteration and imposed the assumption of fast randomization on the remaining trajectories.
As described above, Davis' observation is due a property of the F 1 -F 1 tangle-R 1 is mostly occupied by L 0,1 (0) ∩ L 1,0 (1). We quantify this proportion below.
The assumption of fast randomization of the other trajectories and a 50% probability of them reacting is more difficult to support. From the analysis of lobes we know that however intricate the dynamics is, there is no reason for precisely half of the remaining trajectories to leave to reactants and half to products. Instead we find that 123 for small energies, trajectories that spend 2 and more iterations R 0 and R 0 make up a significant part of the tangles (up to half at 0.02350), but their total proportion is very small and only grows slowly with increasing energy. In the interval up to 0.03000, these trajectories make up at most 3% of the total, 2% below 0.02650, see Table 1. Consequently, any estimate of the reaction rate that takes trajectories escaping after 1 iteration into account is accurate to within 3% below 0.03000 and when we include trajectories escaping after 2 iterations, this number drops to less than 1%.
The difficulty lies in accurately calculating the amount of trajectories. At the cost of accuracy, Davis used VTST as a measure of trajectories entering the interaction region, μ(L 3,0 (0)) to estimate the size of the tangle and μ(L 3,0 (0) ∩ L 0,2 (1)) to subtract trajectories escaping after 1 iteration. The upper and lower estimates assume all, respectively none, of the trajectories that escape after 2 or more iterations are reactive.

The intricate energy interval
The energy interval 0.02215 < E < 0.02654, when TST is not exact and F 0 is a TS, has been largely avoided in the past. The interaction of invariant manifolds of two TSs posed enough difficulties. Dividing tangles using pieces of invariant manifolds and following pips to understand dynamics within make this task possible. We divide tangles differently to the lobe dynamics approach, because we aim to describe and measure parts of heteroclinic tangles that do not necessarily fall into a single lobe.

Division of a tangle
Davis [6] calculated pieces of invariant manifolds in this interval at an energy of 0.7 eV ≈ 0.02572, but complexity of their intersections did not admit deeper insight. With current understanding it is not possible to consider all the invariant manifolds at once, because even identifying lobes is challenging, not to speak of their intersections.
We use the approach outlined in Sect. 5.5 and concentrate on W F 1 and W F 1 , while keeping W F 0 in mind near F 0 . A similar approach may be used for homoclinic tangles. We separate predictably evolving trajectories from chaotic ones, for example trajectories escaping after 1, 2 or 3 iterations from the rest of the tangle. To our knowledge, tools for identifying particular lobe intersections and determining the area, a heteroclinic tangle surgery toolbox, have not been previously presented or reported.
There is one more important property of the manifolds that stands out from all previous figures. Inside the F 1 -F 1 tangle, W u F 1 − and W u F 1 + are restricted to the stripe between two pieces of unstable manifold, e.g. U [ F 1 , Q 3 ] and U [F 1 , Q 1 ] at 0.02700 in Fig. 24 or U [ F 1 , P Q 1 ] and U [F 1 , P Q 1 ] at 0.02400 in Fig. 25. Similarly W s F 1 − and W s F 1 + are confined to a single stripe. We remark that W F 0 are located between W F 1 − and W F 1 + and thereby confined as well. It therefore makes sense to study this stripe in detail.  Fig. 25, and the lower part is symmetric to it. Each lobe consists of two disjoint sets, for example, We remark that lobes do not intersect outside R 5 and leave the interaction region. Disjoint capture lobes imply: Remark 3 R 5 has the no-return property.
As found in Sect. 5.6, a large part of R 5 behaves regularly and leaves the tangle within 1 iteration. As argued in Sect. 5.5, stable manifolds contract in forward time and thereby act as a barrier. Everything above W s F 1 − leaves at the next iteration to reactants, everything below W s F 1 + leaves to products. This agrees with the lobes L 5,3 (1) and L 5,2 (1) that leave R 5 by definition.
The remainder of R 5 is the stripe between W s F 1 − and W s F 1 + , the only part of R 5 where stable manifolds can lie. We refer to it as the capture stripe and denote it R 6 , see Fig. 25. Its boundary consists of S[ F 1 , In backward time, the roles of stable and unstable manifolds switch-everything below W u F 1 − and above W u F 1 + escapes R 5 . Define R 7 , the escape stripe bounded by We conclude that all complicated and chaotic dynamics is confined to R 6 ∩ R 7 and due to the no-return property of R 5 : Remark 4 R 6 and R 7 have the no-return property. Note that the boundary of R 7 is the image of the boundary of R 6 . Necessarily and due to preservation of area μ(R 6 ) = μ(R 7 ).
There are more regions with the no-return property in the F 1 -F 1 tangle. Obviously, R 5 \(R 6 ∪ R 7 ) must be a no-return region as it escapes the R 5 immediately after entering. Also R 6 \R 5 as the entry point to R 7 must have the no-return property as well as capture and escape lobes.

Dynamical properties
To shorten and facilitate the description of reactive and dynamical properties of R 5 , R 6 and R 7 , we introduce the following classification of trajectories.

Definition 7
We call the set of trajectories: directly reactive (D R) if they remain in R 2 or R 3 , directly nonreactive (DN ) if they do not enter the interaction region, captured reactive after n iterations (C R n ) if they react after n iterations in R 5 , captured nonreactive after n iterations (C N n ) if they return to the region of origin after n iterations in R 5 .
Clearly D R and DN never enter R 5 . Following Sects. 5.6 and 6.1, R 5 \(R 6 ∪ R 7 ) is the region of C N 1 and C R 1 is always empty. C R 2 and C N 2 are pass through R 6 \R 7 and R 7 \R 6 and therefore never enter R 6 ∩ R 7 .
This leaves the complicated evolution and chaotic behaviour restricted to R 6 ∩ R 7 . Below 0.02500, R 6 ∩ R 7 consists of 5 squares near F 0 , F 1 , F 1 , F 2 and F 2 . As F 2 and F 2 approach the bifurcation with F 0 , the three squares near them merge into one around 0.02523 when F 21 bifurcates, see Figs. 26 and 27.
Trajectories enter R 5 via R 6 \R 7 and escape via R 7 \R 6 , hence every trajectory crosses R 6 \R 7 and R 7 \R 6 at most once. The same is true for R 5 \(R 6 ∪ R 7 ) consisting of C N 1 . Therefore of R 5 only the size R 6 ∩ R 7 does not reflect the number of trajectories it contains. It follows that the area of R 6 \R 7 , R 7 \R 6 and R 5 \(R 6 ∪ R 7 ) on the surface of section Σ 0 is the same of their images on DS 1 and DS 1 . Figure 28 shows a more detailed partitioning of R 6 and R 7 . Essentially, R 6 is divided into finer stripes by pieces of W s F 1 − and W s F 1 + that are nearly parallel to the boundary. The boundary of R 6 illustrates how the content of the stripe is deformed when mapped into R 7 . It is compressed along the stable manifolds towards the fixed points, e.g.
and stretched along the unstable manifolds away from the fixed points. We remark that the whole highlighted set in R 5 \R 7 of Fig. 28 is connected, only separated by stable invariant manifolds. When mapped forward it is stretched, but remains connected. The sets labeled by yellow, red and green are alternately mapped above and below the capture stripe.
There is a connection between these coloured stripes and lobe intersections, but lobes do not distinguish how often trajectories cross DS 0 , which is necessary to understand overestimation of the reaction rate by TST.
The connected components of R 6 ∩ R 7 contain dynamics similar to Smale's horseshoe dynamics, [12]. As a consequence we observe a fractal structure, as can be seen in Fig. 29. R 6 ∩ R 7 accounts for less than 12% of the all trajectories that pass through R 5 below 0.02400 when the dynamics is relatively slow. The proportion drops to roughly 7% of R 5 at 0.03000 and remains below 1% of the total amount of trajectories.

Areas
In this system, determining the area of R 5 , R 6 , R 7 , R 5 \(R 6 ∪ R 7 ), R 6 \R 7 and R 6 ∩ R 7 , is significantly easier than calculating lobe intersections. We employ a Monte Carlo based method that is expensive, yet simple. Ultimately the cost and accuracy depend on the level of detail in R 6 ∩ R 7 , i.e. it can be determined a priori. We also tune initial and terminal conditions to obtain a high accuracy at a reasonable cost.
Previous works seem to consider initial conditions on r 1 + r 2 2 = 50, p r 1 < 0, which is a surface near q 2 = 1181. We prefer to sample the hemisphere of DS 1 , through which trajectories enter the interaction region. Directly we have that the difference 123 Fig. 28 Coloured sets in R 6 showing how part of the capture stripe is mapped at 0.02400. W F 1 − are drawn with solid lines, W F 1 + are dashed. C N 1 are shown in orange, C R 2 green and yellow, C N 2 red. Part of blue also belongs to C N 2 . Blue, yellow, red and green are separated by white stripes that are mapped to R 6 ∩ R 7 (Color figure online) between the inward hemisphere of DS 1 and r 1 + r 2 2 = 50, p r 1 < 0 corresponds to DN trajectories.
The slowest of D R are located near the boundary of R 5 , and those near pips evolve similarly to pips. Using pips on W u F 1 we define checkpoints, that mark distance these pips are mapped, i.e. the least distance D R cover in the interaction region in 1 iteration. Then all trajectories that pass the second checkpoint 1 iteration after they pass the first checkpoint, are D R and are not captured in R 5 . Recall that all captured reactive trajectories spend at least 2 iterations in R 5 . Trajectories that have a delay of n iteration between crossing of the checkpoints are C R n . Since R 5 is symmetric, we use the symmetric counterpart of the second checkpoint to identify C N n . At 0.02400, Q 1 and P Q 1 are the natural choice for checkpoints, because mark the endpoints of R 6 \R 7 via which trajectories enter R 5 , see Fig. 30. We define checkpoint Ch Q 1 as a vertical line passing through Q 1 . It is necessary that Ch Q 1 avoids capture lobes, therefore at energies above 0.02900 the computationally most efficient solution is to use another vertical line between Q 1 and F 0 .
The role of the second checkpoint, Ch P Q 1 , is to distinguish trajectories in R 5 from those outside R 5 . We use a linear approximation of S[ F 1 , P Q 1 ], the boundary between the escape lobe and R 5 , in conjunction with a vertical line passing through P Q 1 .
The checkpoint symmetric to Ch P Q 1 is defined analogously and denoted Ch P Q 1 . If desired, we can track the number crossings of DS 0 using the sign of q 2 .
We can measure individual components of R 5 : μ(R 5 \R 7 ) corresponds to the number of captured trajectories, μ(R 5 \(R 6 ∪ R 7 )) is given by C N 1 . Then and R 6 ∩ R 7 can be deduced from C R n and C N n where n ≥ 3. The latter follows from the fact that C R 2 and C N 2 do not pass through R 6 ∩ R 7 .
This method is not computationally cheap, but the computational difficulty can be easily estimated a priori. Determining the distribution up to C R n and C N n with N initial conditions requires approximately n N iterations of the map P, but considering the prevalence of D R and DN , this number will be considerably lower.
Alternative approaches to calculating lobe areas face the obstacle in distinguishing the inside from the outside of a lobe, not to mention their intersections. Recent developments [17,18] suggest that a reactive island approach can be used to calculate areas of intersections in a cheaper and simpler manner.

Quantification
In Table 1 we present proportions of areas of classes of trajectories on the plane r 1 + r 2 2 = 50, p r 1 < 0. Between 10 7 and 2.10 8 initial conditions were used to obtain these values.
The proportion of DN decreases steadily over the whole interval presented in Table 1 and beyond. This is not surprising given that widening bottlenecks allow more trajectories enter the interaction region. Thereby nonreactive heavily oscillating trajectories enter the interaction region and consequently R 5 \(R 6 ∪ R 7 ) grows faster than the rest of R 5 .
Note that the proportion of D R culminates between 0.02500 and 0.02550. At this energies the geometry of the F 1 − F 1 tangle simplifies with the consequence that all captured trajectories cross DS 1 . The proportion of D R above 0.02550 decreases predominantly in favour of C N 1 . We observe the growth of capture lobes mainly in the area of large | p 2 | momentum, containing predominantly C N 1 trajectories (Table 1), approaching the maximal values of | p 2 | at the given energy. This implies that for a given (small) | p 1 | momentum, trajectories are more likely to react at a lower energy due to smaller capture lobes.
In the physical world large values of | p 2 | correspond by definition (Sect. 3.3) to large | p r 1 | on the reactant side and large | p r 2 | on the product side. Since trajectories are less likely to react at higher energies, the mechanism for transfer of kinetic energy between the degrees of freedom in the interaction region must be failing at high energies. Consequently, the energy passed from the incoming H to the H 2 may be so high, that it repels the whole molecule instead of breaking its bond. This may be true Directly reactive (DR) and directly nonreactive (DN ) trajectories do not enter R 5 . Captured reactive (C R 2 ) and captured nonreactive (C N 1 , C N 2 ) enter and leave R 5 after 1 or 2 iterations. Other trajectories do not leave R 5 within 2 iterations after their entry and are inside R 6 ∩ R 7 . Horizontal lines represent the creation of the homoclinic tangles and loss of normal hyperbolicity of F 0 for a whole class of collinear atom-diatom reactions, provided it is possible to define an interaction region multiple TSs.

MC based bounds
Using Table 1 we are able to formulate estimates of the reaction rate up to arbitrary precision. The idea is similar to [6]. An upper/lower bound on the reaction rate is obtained by assuming that all/none of the trajectories that remain in R 5 after n iterations react. Table 2 contains the resulting bounds. Denote U 1 the rate estimate obtained by assuming all trajectories in R 5 \(R 6 ∪ R 7 ) react, or equivalently only C N 1 do not react. Since C N 2 and some of Other do not react, the true reaction rate is lower. Lemma 4 U 1 = μ(D R) + μ(C R 2 ) + μ(C N 2 ) + μ(Other) is an upper bound of the reaction rate. U 1 can be easily improved by acknowledging that C N 2 are nonreactive. Denote this bound by U 2 . Because the R 6 ∩ R 7 contain reactive as well as nonreactive trajectories, the true reaction rate is lower. A lower bound L 2 is obtained assuming all of R 6 ∩ R 7 are nonreactive trajectories.

Lemma 6 L 2 = μ(D R) + μ(C R2) is a lower bound of the reaction rate.
The difference between L 2 and U 2 is precisely μ(Other). This gives us an upper bound on the error of both estimates. An estimate of the reaction rate can be obtained using L 2 and U 2 .

Conclusion
We have studied invariant manifolds of TSs to find an explanation for the decrease of the reaction rate. In the process of understanding how energy surface volume passes through homoclinic and heteroclinic tangles formed by these invariant manifolds we found the need for tools that would allow us to work with the tangles and not get lost in details of its chaotic structure. We introduced a suitable division of homoclinic and heteroclinic tangles that is simple and understandable based on reactive properties of trajectories.
Once divided, the heteroclinic tangles decompose into areas of simple and more complicated dynamics. We were able to identify a large class of trajectories that are merely diverted by the tangles and areas of fractal horseshoe-like structure near hyperbolic or inverse-hyperbolic periodic orbits. In addition to a better understanding, the division provides an easy way calculating the corresponding areas.
Contrary to expectations, the decline of the reaction rate is not a result of loss of normal hyperbolicity. We may consider the decrease of the reaction rate and loss of normal hyperbolicity to be consequences of insufficient transfer of kinetic energy between the degrees of freedom. In physical terms, the single atom has so much kinetic energy, that it repels the whole molecule instead of becoming part of it.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.