We examine the phase space structures that govern reaction dynamics in the absence of critical points on the potential energy surface. We show that in the vicinity of hyperbolic invariant tori, it is possible to define phase space dividing surfaces that are analogous to the dividing surfaces governing transition from reactants to products near a critical point of the potential energy surface. We investigate the problem of capture of an atom by a diatomic molecule and show that a normally hyperbolic invariant manifold exists at large atom-diatom distances, away from any critical points on the potential. This normally hyperbolic invariant manifold is the anchor for the construction of a dividing surface in phase space, which defines the outer or loose transition state governing capture dynamics. We present an algorithm for sampling an approximate capture dividing surface, and apply our methods to the recombination of the ozone molecule. We treat both 2 and 3 degrees of freedom models with zero total angular momentum. We have located the normally hyperbolic invariant manifold from which the orbiting (outer) transition state is constructed. This forms the basis for our analysis of trajectories for ozone in general, but with particular emphasis on the roaming trajectories.

## I. INTRODUCTION

Critical points of the potential energy surface (PES) have played, and continue to play, a significant role in how one thinks about transformations of physical systems.^{1,2} The term “transformation” may refer to chemical reactions such as isomerizations^{3–14} or the analogue of phase transitions for finite size systems.^{2,15,16} A comprehensive description of this so-called “energy landscape paradigm” is given in Ref. 2. The energy landscape approach is an attempt to understand dynamics in the context of the geometrical features of the potential energy surface, i.e., a configuration space approach. However, the arena for dynamics is phase space,^{17–19} and numerous studies of nonlinear dynamical systems have taught us that the rich variety of dynamical behavior possible in nonlinear systems *cannot* be inferred from geometrical properties of the potential energy surface alone. An instructive example is the fact that the well-studied and nonintegrable Hénon-Heiles potential can be obtained by series expansion of the completely integrable Toda system.^{20} Nevertheless, the configuration space based landscape paradigm is physically very compelling, and there has been a great deal of recent work describing *phase space signatures* of index one saddles^{21} of the potential energy surface that are relevant to reaction dynamics (see, for example, Refs. 22–25). More recently, index two^{26–28} and higher index^{29} saddles have been studied (see also Refs. 30–32).

The work on index one saddles has shown that, in *phase space*, the role of the saddle *point* is played by an *invariant manifold* of saddle stability type, a so-called normally hyperbolic invariant manifold or NHIM.^{33,34} The NHIM proves to be the anchor for the construction of dividing surfaces (DSs) that have the properties of no (local) recrossing of trajectories and minimal (directional) flux.^{35} There is even richer variety of phase space structures and invariant manifolds associated with index two saddles of the potential energy surface, and their implications for reaction dynamics have begun to be explored.^{27,31,32} Fundamental theorems assure the existence of these phase space structures and invariant manifolds for a range of energy above that of the saddle.^{34} However, the precise extent of this range, as well as the nature and consequences of any bifurcations of the phase space structures and invariant manifolds that might occur as energy is increased, is not known and is a topic of continuing research.^{36–39}

While work relating phase space structures and invariant manifolds to saddle points on the potential energy surface has provided new insights and techniques for studying reaction dynamics,^{22–25} it by no means exhausts all of the rich possibilities of dynamical phenomena associated with reactions. In fact, a body of work has called into question the utility of concepts such as the reaction path and/or transition state (TS).^{40–49} Of particular interest for the present work is the recognition that there are important classes of chemical reaction, such as ion-molecule reactions and association reactions in barrierless systems, for which the transition state is not necessarily directly associated with the presence of a saddle point on the potential energy surface; such transition states might be generated dynamically, and so are associated with critical points of the amended or effective potential, which includes centrifugal contributions to the energy.^{50–52} The phenomenon of transition state switching in ion-molecule reactions^{53–55} provides a good example of the dynamical complexity possible in such systems, as does the widely studied phenomenon of roaming.^{41,46,49}

The lack of an appropriate critical point on the potential energy surface with which to associate a dividing surface separating reactants from products in such systems does *not* however mean that there are no relevant geometric structures and invariant manifolds in *phase space*. In this paper, we discuss the existence of NHIMs, along with their stable and unstable manifolds and associated dividing surfaces, in regions of phase space that do not correspond to saddle points of the potential energy surface. We describe a theoretical framework for describing and computing such NHIMs. Like the methods associated with index one and two saddles, the method we develop for realizing the existence of NHIMs is based on normal form theory; however, rather than normal form theory for saddle-type equilibrium points of Hamilton’s equations (which are the phase space manifestation of index one and two saddles of the potential energy surface), we use normal form theory for certain hyperbolic invariant tori. The hyperbolic invariant tori (and their stable and unstable manifolds) alone are not adequate, in terms of their dimension, for constructing NHIMs that have codimension one stable and unstable manifolds (in a fixed energy surface). However, by analogy with the use of index one saddles to infer the existence of NHIMs (together with their stable and unstable manifolds, and other dividing surfaces having appropriate dimensions), these particular hyperbolic invariant tori can likewise be used to infer the existence of phase space structures that *are* appropriate for describing reaction dynamics in situations where there is no critical point of the potential energy surface in the relevant region of configuration space.

This paper is organized as follows: In Sec. II, we introduce a model for atom-diatom collision dynamics; specifically, we study capture of an atom by passage through an outer or “loose” transition state (OTS). This dividing surface is not associated with a critical point on the potential surface; rather, the relevant phase space object is associated with a NHIM that exists in the vicinity of partially hyperbolic tori. The general mathematical framework for describing phase space structure in the vicinity of partially hyperbolic tori is reviewed in Appendix A. In Sec. II, we describe the associated NHIM and corresponding (codimension one) DS on the energy shell. By exploiting the adiabatic separation of the (fast) diatomic vibration in the phase space region of interest, we are able to formulate a practical algorithm for sampling the relevant phase space DS so that trajectory initial conditions on the DS can be obtained in a systematic way. Section III applies this algorithm to sample the DS in a model for the initial recombination reaction in ozone, O + O_{2} → O_{3}. Both 2 and 3 degrees of freedom (DoF) models are considered. An interesting finding is that roaming dynamics (as defined and explored in our earlier work ^{56–59}) is important in this reaction, and may therefore be of significance in understanding the anomalous isotope effect.^{60–64} Section IV concludes.

## II. A 3 DOF SYSTEM: A 2 DOF SUBSYSTEM WEAKLY COUPLED TO AN ELLIPTIC DOF APPLICABLE TO CAPTURE THEORY IN ATOM-DIATOM REACTIONS

We consider a 2 DoF subsystem weakly coupled to an elliptic DoF (i.e., a 1 DoF nonlinear oscillator). This situation arises, for example, in triatomic molecules when one of the atoms is far from the diatomic fragment and serves as a model for capture of an atom by a diatomic molecule. Specifically, this is the dynamical model for our trajectory study of the ozone recombination reaction (see Sec. III). In Jacobi coordinates, the diatomic molecule vibration adiabatically decouples^{20} from the DoF describing the distance of the third atom to the center of mass of the diatomic and the angle between the diatomic and the third atom.^{65}

In the following, we will first specify this model more precisely and show how a NHIM can exist in the model. In Appendix A, following Ref. 66, we show that, near a partially hyperbolic invariant torus, a NHIM exists from which a codimension one DS can be constructed having the no-recrossing property.^{35} Such an object is relevant for constructing a DS for incoming trajectories in our model.

We also describe a method for sampling the DS for a 2 DoF system and a special 3 DoF system which is relevant to the roaming scenarios that we study.

### A. The existence of a NHIM for a 2 DoF subsystem weakly coupled to an elliptic oscillator

In order to show the existence of a NHIM, we first describe the model to be analyzed. We consider a Hamiltonian describing a triatomic molecule in Jacobi coordinates for total angular momentum *J* = 0 (rotating in the plane) of the following form:

where *r* is the diatomic internuclear distance, *R* is the distance from the centre of mass of the diatomic to the third atom, and *θ* is the angle between *r* and *R*, *p _{r}*,

*p*and

_{R}*p*

_{θ}are the corresponding conjugate momenta, and

*μ*and

_{r}*μ*are reduced masses. At large

_{R}*R*(ratio

*r*/

*R*≪ 1), the (fast) diatomic vibrational mode becomes adiabatically separable from the (

*r*,

*θ*) DoF.

^{20,65}The approximate decoupling of the (

*R*,

*p*) and (

_{R}*θ*,

*p*

_{θ}) DoF can be expressed by assuming that for large atom-diatom distances (small parameter

*r*/

*R*) the potential

*V*(

*r*,

*R*,

*θ*) has the form

For example, using a multipole expansion,^{67} individual terms in the above equation can be written explicitly as

In Eq. (3a), *V _{r}*(

*r*) has been written as a Morse potential representing the interaction between the two atoms forming the diatomic. Setting the exponent

*ν*= 4 in Eq. (3b) corresponds to a charge/induced-dipole interaction, for example. The coupling term Eq. (3c) is a power series in the ratio

*r*/

*R*, which as noted is assumed to be small (large

*R*,

*r*/

*R*≪ 1).

To show the existence of a NHIM, we first consider the completely decoupled case, *V*_{coupling}(*r*, *R*, *θ*) = 0 (the potential is isotropic for *R* → ∞). Hamilton’s equations then have the form

Eq. (4f) implies that *p*_{θ} is a constant of the motion. With *p*_{θ} constant, Eqs. (4a), (4d), (4b), and (4e) represent two uncoupled 1 DoF systems for which the solutions (*r*(*t*), *p _{r}*(

*t*)) and (

*R*(

*t*),

*p*(

_{R}*t*)) can be found separately. Inserting these solutions into Eq. (4c), the solution

*θ*(

*t*) can then be determined. The equations of motion for the uncoupled case (i.e.,

*V*

_{coupling}(

*r*,

*R*,

*θ*) = 0) therefore have the form of 3 uncoupled 1 DoF systems. The DoF (

*θ*,

*p*

_{θ}) is an angle-action pair, describing either oscillation or free rotation.

We assume that in the 2 DoF subsystem corresponding to the coordinates (*R*, *θ*, *p _{R}*,

*p*

_{θ}), we have located an unstable periodic orbit. This unstable periodic orbit is the cartesian product of a hyperbolic fixed point in the (

*R*,

*p*) coordinates, denoted $ ( R \u0304 , p \u0304 R ) $ (the “centrifugal barrier” in the 1 DoF (

_{R}*R*,

*p*) subsystem) and a periodic orbit in the 1 DoF oscillator (

_{R}*θ*,

*p*

_{θ}) subsystem which we denote by (

*θ*(

*t*),

*p*

_{θ}= constant). From the discussion in Appendix A, we see that this periodic orbit is a NHIM for the 2 DoF subsystem.

The NHIM for the 3 DoF uncoupled system for a fixed energy *E* is obtained by taking the cartesian product of the unstable periodic orbit in the 2 DoF subsystem with a periodic orbit in the (*r*, *p _{r}*) subsystem, which we denote by ($ r \u0304 ( t ) $, $ p \u0304 r ( t ) $). More precisely,

where *B* ⊂ ℝ denotes the domain of *p*_{θ}. The NHIM is 3-dimensional in the 5-dimensional energy surface. More details of the structure of the NHIM are given in Section II B.

We can now appeal to the general persistence theory for NHIMs under perturbation^{34} to conclude that the NHIM just constructed persists when the term *V*_{coupling}(*r*, *R*, *θ*) is non-zero. This completes the proof of the existence of a NHIM for a diatomic molecule weakly interacting with an atom for total angular momentum *J* = 0.

#### 1. Decoupling in atom-diatom systems: Adiabatic separability

The mathematical results we use concerning NHIMs, etc., which are applied to construct a (partial) DS, require *only* that the diatomic vibration be approximately separable at large *R* (in the vicinity of the centrifugal barrier). The ultimate origin of the separability assumed at large *R* (that is, for ratio *r*/*R* ≪ 1) is to be found in a separation of timescales: the diatomic vibration frequency is much larger than other dynamical frequencies in the relevant region of phase space, so that in the region of phase space studied here the diatomic vibration can be adiabatically separated from the other degrees of freedom.^{20,65,68}

In the context of our study, the breakdown of the relevant adiabatic approximation^{68} is unlikely to occur in the region of phase space of interest. This is because, in order to reduce the vibrational frequency of the diatom substantially, it is necessary to provide sufficient energy to bring the diatomic fragment near dissociation. If we require that the other bond remain extended (i.e., *R* ≫ 1), then such a configuration is near the threshold for collision-induced dissociation, which is well above the energies studied in our paper. In any case, transfer of large amounts of energy from the *R* DoF to the *r* DoF requires that the system pass through either of the tight TSs; the details of the resulting reaction dynamics are beyond the scope of the present study.

The fact that the diatomic vibration does decouple at large *R* means that we can use the NHIM to define a DS which is used to sample trajectory initial conditions. These trajectories are then integrated numerically and it is found, in the case of ozone (see Sec. III), that there is a significant degree of decoupling over the whole phase space region of interest (bounded by the OTS and the 2 inner tight TS, roughly speaking between *R* = 2 Å and *R* = 9 Å). We are therefore able to establish a connection between 2 and 3 DoF “roaming” dynamics (see Sec. III). This latter separability is nevertheless by no means required in order for us to construct and sample the DS associated with the OTS.

### B. Approximate numerical procedure for sampling the DS

In Subsection II A, we demonstrated the existence of a NHIM for a diatomic molecule weakly interacting with an atom. The assumptions used to demonstrate the existence of the NHIM are however not sufficient to provide a straightforward procedure to sample the DS attached to this NHIM. We now make an even stronger assumption that will enable us to sample the DS, albeit approximately.

We start from the fact that for sufficiently large values of *R*, the diatomic vibration is completely decoupled from the other DoF. This means that we assume that (for large *R*, *r*/*R* ≪ 1) the coupling term *V*_{coupling}(*r*, *R*, *θ*) in the potential function Eq. (2) will depend only on *R* and *θ*, i.e., *V*_{coupling}(*r*, *R*, *θ*) → *V*_{coupling}(*R*, *θ*). Again, this decoupling results from an adiabatic separation of the diatomic vibrational mode from the other DoF; the effective potential then depends only on the value of vibrational action, assumed to be constant.

For simplicity, in the 2 DoF subsystem Hamiltonian, we can consider the coordinate *r* to be fixed at its equilibrium value *r _{e}*, i.e., at the minimum of the Morse like potential well. (Alternatively, we can fix the value of

*r*at its average value 〈

*r*〉 associated with diatomic vibration at fixed action.) Under these stronger assumptions, our system becomes a completely decoupled 2 DoF subsystem plus an elliptic oscillator. The sampling of the DS then reduces to a sampling of the part of the DS for the 2 DoF subsystem and a sampling of the other part corresponding to the elliptic oscillator.

In the following, we first describe an algorithm for sampling a DS attached to a NHIM in a 2 DoF system. Then we use this algorithm to develop a procedure for sampling a DS for a system composed of the cartesian product of a 2 DoF subsystem with an elliptic oscillator.

#### 1. Sampling of a dividing surface attached to a NHIM for a 2 DoF system for a fixed total energy

When the theory developed in Appendix A is applied to a 2 DoF system, the hyperbolic torus reduces to an unstable periodic orbit. In other words, the NHIM in this case is just an unstable periodic orbit. We now describe an algorithm to sample points on a DS attached to a periodic orbit for a 2 DoF system.

For a 2 DoF system, the DS is topologically homeomorphic either to a 2-sphere of which the NHIM (the PO) is an equator or to a 2-torus (see below). To sample points on the DS, we need to specify the values of four phase space coordinates. The PO is generally obtained by a numerical search.^{69} The NHIM is represented by discrete points in phase space on the PO at different times corresponding to the discretization of the period of the PO. Once projected onto configuration space, these points lie on a curve. This curve in configuration space can be approximated by a spline interpolation in order to obtain a continuous representation of the curve. The sampling of the DS starts by sampling points along this curve in configuration space. This fixes the values of the 2 configuration space coordinates. To determine the values of the remaining 2 phase space coordinates (momenta), we use the fact that the sampling is at a fixed total energy. So, for each point on the projection of the PO into configuration space, we can scan the value of one of the momenta and then obtain the remaining momentum by solving numerically the equation *H* = *E*, *H* being the Hamiltonian of the 2 DoF system, for the remaining momentum coordinate. The maximum value of the first momentum coordinate to be scanned can be determined by setting the second momentum coordinate to zero and solving the equation *H* = *E* for each point of the PO in configuration space.

To specify such an algorithm to sample points on the DS, we denote configuration space coordinates by (*q*_{1}, *q*_{2}) and momentum space coordinates by (*p*_{1}, *p*_{2}). (This notation emphasizes that our approach is applicable to systems other than the atom-diatom model considered here.) The algorithm to sample the DS for a 2 DoF system can be summarized as follows:

Locate an unstable PO.

Perform a spline interpolation

^{70}of the curve obtained by projecting the PO into configuration space.Sample points on that curve and obtain points (

*q*_{1i},*q*_{2i}), for*i*= 1, …,*N*, where_{p}*N*is the desired number of points._{p}For each point, (

*q*_{1i},*q*_{2i}) determine*p*_{1max}by solving*H*(*q*_{1i},*q*_{2i},*p*_{1}, 0) =*E*for*p*_{1}.For each point, (

*q*_{1i},*q*_{2i}) sample*p*_{1}from zero to*p*_{1max}and solve the equation*H*(*q*_{1i},*q*_{2i},*p*_{1j},*p*_{2}) =*E*to obtain*p*_{2}.

For kinetic energies that are quadratic in the momenta, the momentum values satisfying the energy condition *H* = *E* above will appear with plus or minus signs. Computing the points of the DS with both positive and negative for *p*_{1} covers both hemispheres of the DS.

For a DS is obtained by this sampling procedure, we now address the following questions:

What is the dimensionality of the sampled surface?

What is the topology of the sampled surface?

Does the NHIM (PO) bound the 2 hemispheres of the DS?

Do the trajectories cross the DS and what does “reaction” mean with respect to the sampled DS?

Consider the Hamiltonian

where the potential is periodic in coordinate *q*_{2} with period 2*π*.

As we have noted, the sampling procedure begins with the location of an unstable PO, where the projection of the located PO is 1-dimensional. We now turn our attention to the questions listed above.

##### a. Dimensionality.

The dimensionality of the sampled DS is relatively straightforward to determine. As described above, we locate an unstable PO and perform a spline interpolation for the set of points obtained by projecting the PO into configuration space. We then sample points along this curve. For each (*q*_{1}, *q*_{2}) point, we determine the maximum possible value of one of the momenta by solving the equation *H*(*q*_{1}, *q*_{2}, *p*_{1}, 0) = *E* for *p*_{1}, for example. We obtain the maximum value for *p*_{1}, *p*_{1max}, permitted by the energy constraint. We then sample the momentum *p*_{1} from −*p*_{1max} to +*p*_{1max}.

Topologically, the space sampled in configuration space is a 1D segment and the space sampled in the momentum is also a 1D segment. The full space sampled is therefore a 2D surface embedded in the 3D energy surface in the 4D phase space.

##### b. Topology of the DS.

We must distinguish 2 cases, which differ in the type of PO sampled. There are 2 types of PO to consider. The first type is a PO for which the range of the PO in the periodic coordinate *q*_{2} is less than the full period of the coordinate. The second type is a PO which makes a full cycle in the periodic coordinate.

###### i. Type 1 PO.

This type of PO has a range in the periodic coordinate that is less than the period of the coordinate. These POs have 2 turning points at which the total energy is purely potential and the momenta are zero.

If, for example, we sample the *q*_{2} coordinate, this also fixes *q*_{1}. On a fixed energy surface, we have to satisfy the equation *H*(*q*_{1}, *q*_{2}, *p*_{1}, *p*_{2}) = *E*, so having fixed *q*_{1} and *q*_{2}, the potential energy is fixed at a particular value *V*_{0}. The equation to be satisfied by the momenta is then

This is the equation of an ellipse in the (*p*_{1}, *p*_{2}) plane. Moreover, the lengths of the large and small axes of the ellipse tend to zero as we approach the turning points where *E* = *V*_{0}. So, as we follow the curve of the PO in configuration space from one turning point to the other, the intersection of the DS with the (*p*_{1}, *p*_{2}) plane begins as a point and then becomes a family of ellipses whose 2 axes depend on the potential energy value, eventually shrinking down to a point again as the other turning point is reached. The topology of the DS is then the product of the 1D segment in configuration space with a family of ellipses degenerating to points at either ends, that is, a 2-sphere.

###### ii. Type 2 PO.

This type of PO makes a full cycle in the periodic coordinate. It is important to note that for this case, the DS is constructed from 2 POs. For a PO that makes a full cycle in the periodic coordinate *q*_{2}, the conjugate momentum *p*_{2} never vanishes and always has the same sign (either positive or negative). For a given PO, there is a twin PO for which the momentum has the opposite sign (the PO trajectory traverses the same path but in the opposite direction).

When we sample the curve in configuration space and analyse the shape of the DS in momentum space for a particular point in configuration space, we solve the same equation as for the previous case. Therefore, the shape of the DS in momentum space is also an ellipse. However, now the ellipses never degenerate to points as there are no turning points along the PO at which the total energy is purely potential. Also, since the coordinate *q*_{2} is periodic, the end point of the PO has to be identified with the starting point. In this situation, the DS is then the cartesian product of a circle with a family of ellipses and has the topology of a 2-torus.

##### c. The NHIM separates the DS into two parts.

To discuss this point we also must distinguish between type 1 and type 2 POs.

###### i. Type 1 PO.

As above, we fix a particular point in configuration space (*q*_{10}, *q*_{20}). As we scan one of the momenta, there will be a particular value corresponding to the momentum value on the PO. Since, the sampling is at the same energy as the PO energy, the value of the remaining momentum will also correspond to the value on the PO. For the configuration space point (*q*_{10}, *q*_{20}), this happens exactly twice. For a nearby configuration space point $ ( q 10 \u2032 , q 20 \u2032 ) $, the 2 momentum space points belonging to the PO will be near the 2 points for the point (*q*_{10}, *q*_{20}).

As we move along the PO in configuration space, these 2 points in momentum space will trace out the PO in phase space. At the turning points, the ellipse in momentum space shrinks to a point which is also the meeting point of the 2 points in momentum space belonging to the PO. So, the NHIM (the PO) belongs to the 2D sphere we sample. It forms a 1D closed curve which divides the sphere into 2 hemispheres.

###### ii. Type 2 PO.

For this type of PO, if we fix a configuration space point (*q*_{10}, *q*_{20}) and scan one of the momentum variables, we will meet a point for which the momentum equals the value of that momentum on the PO. The other momentum has its value equal to its value on the PO due to the energy constraint. We will also encounter a point for which the value of the momentum we are scanning is equal to the value of that momentum coordinate on the twin PO (the one which moves in the opposite direction with *p*_{2} having the opposite sign of the original PO). The corresponding value of the remaining momentum coordinate is also on the twin PO.

The 2 POs belong to the surface we are sampling. Taken individually, these two curves do not divide the 2-torus into two parts but the DS is divided into two parts when we consider both POs together: one part of the DS connects the *p*_{2} > 0 PO with the *p*_{2} < 0 PO on the left of the ellipse in momentum space and one part connects the *p*_{2} > 0 PO with the *p*_{2} < 0 PO on the right of the ellipse, see Fig. 1.

##### d. Trajectories cross the DS and the meaning of reaction for the constructed DS.

Except for points on the NHIM, which belongs to the surface we are sampling, no trajectory initiated on the DS remains on the DS. Off the NHIM, the Hamiltonian vector field is everywhere transverse to the sampled surface, and every trajectory initiated on the DS has to leave the surface, i.e., trajectories cross the DS.

We have shown that the NHIM (or the 2 NHIMs for type 2 POs) divides the DS into two parts. One part of the DS intersects trajectories evolving from reactants to products and the other intersects trajectories evolving from products to reactants. Here, “reaction” means crossing one of the two parts of the DS and the direction of the reaction, forward or backward, depends on which half of the DS the trajectory crosses.

To show that the vector field is nowhere tangent to the sampled surface, except on the NHIM, we must evaluate the flux form^{27} associated with the Hamiltonian vector field through the DS on tangent vectors to the DS. The Hamiltonian vector field is not tangent to the DS if the flux form is nowhere zero on the DS, except at the NHIM. Details of the calculation to verify this assertion are given in Appendix B.

#### 2. Sampling of a dividing surface attached to a NHIM for a 2 DoF subsystem plus an elliptic oscillator at a fixed total energy

Under the assumptions introduced at the beginning of this section, the problem of the capture of an atom by a diatomic molecule at zero total angular momentum reduces to consideration of the dynamics of a 2 DoF subsystem and an elliptic oscillator. The algorithm presented for the sampling of a DS attached to a PO is an essential element in our derivation of a sampling procedure of the DS in the case of a 3 DoF system composed of a 2 DoF subsystem plus an elliptic oscillator.

To begin, we need to understand the nature of the NHIM for this case. The motion for the elliptic oscillator (the diatomic vibration) is periodic. For sufficiently small amplitude vibration (for example, at energies below the dissociation energy of the Morse like potential), these POs are stable. The NHIM for the full 3 DoF system (2 DoF subsystem plus 1 DoF elliptic oscillator) is obtained by considering the direct product of the NHIM for the 2 DoF subsystem with a PO of the elliptic oscillator. The NHIM for the 2 DoF subsystem, as we saw above, is just an unstable PO. If we consider a particular unstable PO for the 2 DoF subsystem and a particular PO for the elliptic oscillator, the structure obtained by the cartesian product of these 2 manifolds is a 2-dimensional torus. Every particular partitioning of the total energy between the 2 DoF subsystem and the elliptic oscillator corresponds to a particular 2-dimensional torus for the NHIM. The NHIM is then made up of a 1-parameter family of 2-dimensional tori, where the parameter of this family determines the distribution of energy between the 2 DoF subsystem and the elliptic oscillator.

Now that we have the structure of the NHIM, it remains to construct the associated DS. Like the NHIM, the DS has a product structure: one part is associated with the 2 DoF subsystem and the other with the elliptic oscillator. For the elliptic oscillator, all phase space points belong to a particular PO. For the 2 DoF subsystem, the component of the DS consists of a DS attached to a PO as we described in Sec. II B. If we consider one particular partitioning of the total energy *E* between the 2 DoF subsystem and the elliptic oscillator, the portion of the DS associated with this distribution consists of the direct product of the DS associated with the unstable PO for the 2 DoF subsystem and the PO for the elliptic oscillator. The full DS is obtained by considering the DS for all possible partitions of the energy. The full DS is then foliated by “leaves,” where each leaf of the DS is associated with a particular partitioning of the energy.

To implement an algorithm for sampling, the DS we need to determine the potential *V _{r}*(

*r*) for the diatomic molecule. This can be done by fixing a sufficiently large value of

*R*, and fixing a value for

*θ*. In general, for large values of

*R*, the potential

*V*(

*r*,

*R*,

*θ*) is isotropic, so any value of

*θ*is sufficient. One can then use spline interpolation to obtain an analytical representation of the resulting Morse-like potential.

An algorithm to sample the DS in the case of a decoupled 2 DoF subsystem plus an elliptic oscillator is given as follows:

Consider a total energy

*E*for the full system.Consider a particular partitioning of the total energy

*E*between the 2 DoF subsystem and the elliptic oscillator.Sample points on the PO corresponding to the energy associated with the elliptic oscillator.

Sample points on the DS attached to the PO corresponding to the energy associated with the 2 DoF subsystem using the algorithm presented in the Sec. II B.

Take the direct product the two sets of points obtained in the two previous steps.

Repeat the previous steps for many different energy distributions between the 2 DoF subsystem and the elliptic oscillator.

There is a subtlety related to the sampling procedure just described. The hyperbolic POs of the 2 DoF subsystem are associated with the centrifugal barrier arising in the atom-diatom Hamiltonian in Jacobi coordinates. Following the family of these POs, it is found that as the energy decreases to the dissociation threshold, their location in the (*R*, *θ*) plane moves to larger and larger values of *R*. This suggests that the family goes asymptotically to a neutrally stable PO located at *R* → ∞.^{71} Reducing the energy in the 2 DoF subsystem therefore means locating POs at larger and larger distances in *R*, with longer and longer periods. Location of all such POs is clearly not possible numerically so that the sampling procedure must stop at some finite (small) energy in the diatomic vibration. The sampling procedure will nevertheless cover the DS, at least partially.

Alternatively, the energy in the diatomic vibration can simply be fixed at a specific value, as in standard quasiclassical trajectory method (see, for example, the work of Bonnet and Rayez^{72} in which the energy of the diatomic vibration is fixed at the zero point energy). This approach then samples a particular “leaf” of the classical DS.

## III. AN EXAMPLE: RECOMBINATION OF THE OZONE MOLECULE

The recombination reaction O + O_{2} → O_{3} in the ozone molecule is known to exhibit an unconventional isotope effect.^{60} There is a large literature on this subject and several theoretical models have been proposed to identify the origin of this isotope effect.^{60,61,64} One key to understanding the isotope effect in ozone is to model the recombination of an oxygen molecule with an oxygen atom. Marcus and co-workers have made a detailed study of the isotope effect in ozone^{61–63} and have used different models for the transition state in the recombination of ozone (tight and loose TS). In their studies, they used a variational version of Rice-Ramsperger-Kassel-Marcus (RRKM) theory^{73,74} (Variational Transition State Theory (VTST)) to locate the transition state. In order to obtain agreement between their calculations and experimental results, they introduced parameters whose role consisted of correcting the statistical assumptions inherent in RRKM theory. Marcus and co-workers proposed several explanations for the physical origin of these parameters^{61–63} but their significance is still not fully understood. The deviation from statistical behavior in the ozone molecule suggests that a dynamical study of the recombination reaction of ozone rather than a statistical model might shed light on the origin of the unconventional isotope effect. Such a study necessitates the use of DSs having dynamical significance, i.e., constructed in phase space.

To illustrate the construction of phase space DS for the capture of an atom by a diatomic molecule, we apply the algorithms described in Sec. II to the ozone recombination reaction. To attack this problem we require a PES for the ozone molecule which describes large amplitude motion. Schinke and co-workers produced the first accurate PES for ozone.^{75,76} This PES exhibits a small barrier in the dissociation channel just below the asymptote for dissociation forming a van der Waals minimum and an associated “reef structure.” Recent *ab initio* calculations^{77–79} have however shown that the barrier is not present, and that the potential increases monotonically along the dissociation coordinate (there is no saddle point along the reaction coordinate). Here, we use the PES produced by Tyuterev and co-workers,^{77} which has recently been used to interpret experimental results.^{80}

In the following, we treat the recombination of ozone molecule for the isotopic combination ^{16}O^{16}O^{16}O using two models. First, we treat the problem with a reduced 2 DoF model with zero total angular momentum, and then we look at the 3 DoF problem also with zero angular momentum.

### A. 2 DoF model

The 2 DoF case is easily treated using the analysis of Sec. II A by freezing the diatomic vibration, i.e., (*r*, *p _{r}*) variables. We therefore fix the diatomic bond length at its equilibrium value,

*r*, and set

_{e}*p*= 0. The dependence on

_{r}*r*in the coupling term of the potential then drops out and we have simply

*V*

_{coupling}(

*R*,

*θ*) =

*V*

_{coupling}(

*r*=

*r*,

_{e}*R*,

*θ*). The NHIM is obtained as previously by considering the uncoupled case with

*V*

_{coupling}(

*R*,

*θ*) = 0 and a hyperbolic equilibrium point in the (

*R*,

*p*) DoF, and is just a periodic orbit in the (

_{R}*θ*,

*p*

_{θ}) DoF. This NHIM persists for non-zero coupling and consists of a periodic orbit where variables (

*R*,

*p*) and (

_{R}*θ*,

*p*

_{θ}) are in general coupled.

The sampling of the DS follows the algorithm presented in Sec. II B 1. We use the sampled points on the DS to initiate trajectories and follow the approach of the atom to the diatomic molecule. The resulting DS is a 2-dimensional surface in phase space whose points belong to the chosen energy surface. For our simulation, we choose the total energy *E* = 9200 cm^{−1}. The DS constructed from our algorithm is the first “portal” the system has to cross in order for the atom to interact with the diatomic molecule. This DS is then the phase space realisation of the so-called OTS.^{50,52} At the so-called Tight Transition State (TTS),^{56} the coordinate *θ* undergoes small oscillations. In our phase space setting, these TTSs for a 2 DoF system are DSs attached to unstable POs. In the following, the terms OTS and TTS are to be understood as synonymous with DSs constructed from NHIMs (POs for a 2 DoF system). An important mechanistic issue arising in modeling a reaction is to determine whether the reaction is mediated by an OTS or a TTS. However, in many situations, both types of TS co-exist for the same energy and one has to understand the dynamical implications of the presence of these two types of TS.

In a recent work, we analysed this question for a model ion-molecule system originally considered by Chesnavich.^{53–55} We established a connection between the dynamics induced by the presence of these 2 types of TS and the phenomenon of roaming.^{56,57} Specifically, our phase space approach enabled us to define unambiguously a region of phase space bounded by the OTS and TTS which we called the *roaming region*. We were also able to classify trajectories initiated on the OTS into four different classes. These four classes were derived from two main categories. The main category consists of reactive trajectories. These trajectories start at the OTS and enter the roaming region between the OTS and TTS. They can then exhibit complicated dynamics in the roaming region and eventually react by crossing the TTS to form a bound triatomic molecule. Within this category, we distinguish between *roaming* trajectories making some oscillations in the roaming region (oscillations were defined by the presence of turning points in *R* coordinate) and *direct* trajectories that do not oscillate in the roaming region. The other main category of trajectories consists of trajectories that do not react and for which the atom finally separates from the diatomic molecule. Again we distinguish between trajectories which “roam” (making oscillations in the roaming region) and those that escape directly by bouncing off the hard wall at small values of *R*.

For the ozone recombination problem reduced to a 2 DoF model, we employ the same classification. Figure 2 shows the four classes of trajectories obtained by propagating incoming trajectories originated at the OTS. The black line at *R* ≃ 8 Å is the PO from which the OTS is constructed. The two black curves at *R* ≃ 2.5 Å are the two POs from which the two TTS are constructed. These two TTSs correspond to the two channels by which the ozone molecule can recombine. For the case of ^{16}O_{3}, the two channels reform the same ozone molecule but for different isotope composition, these channels may lead to different ozone molecules. These trajectories show that the roaming phenomenon is at play in the ozone recombination problem. It would be very interesting to investigate how the roaming dynamics varies as the masses of the different oxygen atoms change, as this may provide new insight into the isotope effect in ozone recombination.

### B. 3 DoF Model

We now consider the problem of ozone recombination in a 3 DoF model with zero total angular momentum. Addition of angular momentum DoF is in principle possible in our phase space description of reaction dynamics, but consideration of these DoF will increase the dimensionality of the problem and will require a more elaborate sampling procedure. Consideration of angular momentum in capture problems has been investigated by Wiesenfeld and co-workers^{50–52} and recently by MacKay and Strub.^{81}

The construction of the DS for the OTS follows the procedure described in Sec. II B 2. As explained in Sec. II B 2, we cannot sample the full DS. In our simulation, we sample one leaf of the DS by considering a single partitioning of the energy between the two DoF subsystem and the diatomic vibration. Here, for the sampling of the DS, we fixed the diatomic length *r* to its equilibrium value and use the same value of the total energy as for the 2 DoF model, 9200 cm^{−1}.

The sampled points on the DS were again used to propagate trajectories, which are classified as previously into four classes. The result of the propagation is shown on Fig. 3. Again the black line at *R* ≃ 8 Å represents the projection of the PO defining the OTS, while the two black curves at *R* ≃ 2.5 Å are projections of POs. These POs are not strictly (projections of) the proper phase space structures from which TTSs can be constructed, as for a 3 DoF system the NHIM supporting a DS consists of a 1-parameter family of invariant 2-tori. Nevertheless, these POs give an idea of the approximate location of the projection of these structures into configuration space for ozone recombination. For the ozone problem, it nevertheless appears that they are very good approximation to the true NHIMs for the TTSs as we can see in Fig. 3 that none of the non reactive trajectories (panels (c) and (d)) extends beyond the projections of those POs on configuration space. The reason for that is presumably that for ozone the diatomic vibration is well decoupled from the 2 other DoF throughout the entire roaming region. This fact validates the reduced 2 DoF model presented above.

### C. Role of separatrices in the 2D and 3D model

The approach taken in the present paper has been to classify the behavior of trajectories initiated on the OTS, especially with respect to “roaming” behavior. Roaming trajectories are trapped as a consequence of energy transfer between radial (*R*) and angular (*θ*) DoF. This energy transfer could (at least in 2D) be described and analyzed in terms of passage across a broken separatrix following the work of Davis, Gray, and Skodje (see, for example, Refs. 5 and 82). This kind of analysis has already been attempted: see, for example, the paper of Ramachandran and Ezra.^{83} The latter work suggests that the separatrix/turnstile view, while formally applicable, is not necessarily very useful or practical in the context of problems with (internal) rotational degrees of freedom, for the reason that the return time to the surface of section by means of which the separatrix is visualized can diverge; this makes the lobes and turnstiles difficult to plot and visualize.

In the present paper, we have chosen to investigate the dynamics of trajectories entering a region of phase space bounded by well-defined bottlenecks/DS and to classify them with respect to reactive/non-reactive, roaming/non-roaming behavior, as in our other recent studies.^{56–59} The separatrix perspective is complementary to this approach, but more detailed comparison will require further work.

It should nevertheless be noted that in a recent paper,^{59} we have used the lobe/turnstile construction to provide a phase space rationale for the roaming behavior seen in a model for formaldehyde dissociation.

## IV. CONCLUSION

In this paper, we have examined the phase space structures that govern reaction dynamics in the absence of critical points on the PES. We showed that in the vicinity of hyperbolic invariant tori, it is possible to define phase space dividing surfaces that are analogous to the dividing surfaces governing transition from reactants to products near a critical point of the PES.

We investigated the problem of capture of an atom by a diatomic molecule and showed that a NHIM exists at large atom-diatom distances, away from any critical points on the potential. This NHIM is the anchor for the construction of a DS in phase space, which is the entry “portal” through which the atom has to pass in order to interact with the diatomic molecule. This DS defines OTS governing capture dynamics.

Exploiting adiabatic separability of the diatomic vibrational mode in the phase space region of interest, we presented an algorithm for sampling an approximate capture DS. As an illustration of our methods, we applied the algorithm to the recombination of the ozone molecule. We treated both 2 and 3 DoF models with zero total angular momentum. The co-existence of the OTS and TTS in ozone recombination means that roaming dynamics is observed for this reaction. Such roaming dynamics may have important consequences for the unconventional isotope effect in ozone formation.

## Acknowledgments

F.M., P.C., and S.W. acknowledge the support of the Office of Naval Research (Grant No. N00014-01-1-0769), the Leverhulme Trust. B.K.C., F.M., P.C., and S.W. acknowledge the support of Engineering and Physical Sciences Research Council (Grant No. EP/K000489/1). G.S.E. and Z.C.K. acknowledge the support of the National Science Foundation under Grant No. CHE-1223754. We thank Professor Vladimir Tyuterev for providing us with his PES for ozone.

### APPENDIX A: REACTION DYNAMICS IN THE ABSENCE OF EQUILIBRIUM POINTS

Here, we discuss the work of Ref. 66 on partially hyperbolic tori in Hamiltonian systems and associated phase space structures. The starting setup is an *m* DoF canonical Hamiltonian system defined on a 2*m*-dimensional phase space, which we can take to be ℝ^{2m}. This discussion provides the theoretical foundation for our analysis of the atom-diatom system given in Sec. II.

The steps involved in describing phase space structure in the vicinity of partially hyperbolic tori are as follows:

*Step 1: Locate a partially hyperbolic invariant torus in the original system.* First, we define the term “partially hyperbolic invariant torus.” A hyperbolic periodic orbit provides the simplest example. The mathematics literature contains rather technical discussions of the notion of partial hyperbolicity,^{66} but a practical definition is that, under the linearized dynamics, in directions transverse to the torus *and* the directions corresponding to the variables canonically conjugate to the angles defining the torus (this is where the qualifier “partially” comes from), there is an equal number of exponentially growing and exponentially decaying directions. The following steps express the original Hamiltonian in coordinates defined near the torus.

*Step 2: Express the Hamiltonian in coordinates near the partially hyperbolic invariant torus.* Following Ref. 66 we denote the partially hyperbolic *n*-dimensional invariant torus by *N*. In Theorem 5 of Ref. 66 (with the explicit construction of coordinates in Sec. 5), Bolotin and Treschev show that there exist canonical coordinates in a neighborhood of *N* in which the Hamiltonian takes the following form:

where 〈 ⋅ , ⋅ 〉 denotes the usual Euclidean inner product on ℝ^{n}, *A* denotes an *n* × *n* symmetric matrix, and *λ*(*θ*) > 0 for all *θ* ∈ 𝕋^{n} (this is the condition insuring hyperbolicity). The corresponding Hamiltonian vector field is

Note that in these coordinates the invariant torus is given by

Several comments are now in order.

It is important to realize that (A1) is not the “normal form” that is usually discussed in Hamiltonian mechanics.

^{84}This is because*λ*(*θ*) is not necessarily constant. A normal form requires additional transformations to transform*λ*(*θ*) to a constant (the Floquet transformation). We will see that such a transformation is not necessary for constructing NHIMs and DSs.From the form of Equations (A2) the origin of the term “partially” hyperbolic is clear: not all directions transverse to the torus are hyperbolic (the transverse directions that are not hyperbolic are those described by the coordinates

*I*).For

*n*= 1,*N*is a periodic orbit.- For
*n*> 1 the frequency vector of*N*,*ω*, must be*nonresonant*. More precisely, it is required to satisfy a diophantine condition,^{18}i.e.,where $|k|\u2261 \u2211 i = 1 n | k i |$.$ | \u3008 \omega , k \u3009 | > \alpha | k | \u2212 \beta , \alpha , \u2009 \beta > 0 , \u2009 k \u2208 Z n \u2212 { 0} , $ *z*_{−}and*z*_{+}are the exponentially decaying and growing directions, respectively, normal to*N*; they are of equal dimension. Here, we are only considering the case where they are 1-dimensional. This is the analogous case to index one saddles for equilibria. Both*z*_{−}and*z*_{+}can each have dimension greater that one. In that case, we would have a generalization of the notion of index*k*saddles,^{31}for*k*≥ 1, to periodic orbits and tori. We will not consider that case here.

It is useful to emphasize at this point that the phase space is 2*n* + 2 ≡ 2*m* dimensional.

Before using the Hamiltonian (A1) to construct NHIMs and DSs, we (symplectically) rotate the saddle coordinates in the usual way,

to obtain the Hamiltonian,

*Step 3: Construction of phase space structures.* We consider the 2*m* − 1 dimensional energy surface *H*(*θ*, *I*, *q*, *p*) = *h*. The motivation for rotating the saddle in Eq. (A5) was to clarify the meaning of reaction — it corresponds to a change in sign of the *q* coordinate. Therefore, neglecting the $O ( 3 ) $ terms, setting *q* = 0, the forward and backward DSs are given by

respectively, and these two DSs meet at the NHIM, *q* = 0, *p* = 0,

These are interesting equations. They show that the DSs vary with *θ*, but that the NHIM is constant in *θ*, as follows from the form of Hamilton’s equations in these coordinates. At the origin in the space of the saddle degrees of freedom, we have an integrable Hamiltonian system in action-angle variables. We do not attempt to make the hyperbolic DoFs constant by carrying out a Floquet transformation. As a result, the DS (which incorporates hyperbolic DoFs) varies with *θ*. In some sense, one could view the procedure that gives these coordinates as a type of “partial normalization” where we normalize the elliptic DoFs and leave the hyperbolic DoFs alone.

In this appendix, we have shown that near a partially hyperbolic invariant torus a NHIM exists from which a codimension one DS can be constructed having the no-recrossing property.^{35} We have not provided a method for explicitly constructing the coordinates in which we constructed the NHIM and DS. In Section II, we describe a method for sampling the DS for a 2 DoF system and a special 3 DoF system which is relevant to the roaming scenarios that we study.

### APPENDIX B: VERIFICATION THAT TRAJECTORIES ARE TRANSVERSE TO THE DS ASSOCIATED WITH THE NHIM

To carry out this calculation, we require a parametrization of the DS that will enable us derive expressions for tangent vectors to the DS. We can then compute the flux form and evaluate it on tangent vectors to the DS.

*Parametrization of the DS*. The DS is a 2-dimensional surface and therefore we require 2 parameters to parametrize it. These 2 parameters appear naturally from the sampling method, where we sample in the periodic coordinate *q*_{2} and then in the momentum coordinate *p*_{2}. So, we use parameters (*p*_{2}, *q*_{2}). Spline interpolation gives the function *q*_{1} = *q*_{1}(*q*_{2}) describing the dependence of coordinate *q*_{1} on *q*_{2} in the projection of the NHIM onto configuration space. The remaining momentum coordinate, *p*_{1}, is obtained by imposing the constant energy constraint *H* = *E*. The parametrization Φ of the DS can then be specified as

with $1/ I 1 =1/I+1/ ( \mu q 1 2 ) $.

Tangent vectors to the DS are obtained by differentiating the parametrization with respect to the parameters *p*_{2} and *q*_{2},

*The flux form associated with the Hamiltonian vector field*. For a 2 DoF system the phase space volume form Ω is expressed using the symplectic 2-form *ω* as follows:

The energy surface volume 3-form *η* is defined by

The Hamiltonian vector field is

The symplectic 2-form *ω* when applied to an arbitrary vector *ξ* and the vector field *X _{H}* gives

The flux 2-form *φ* through a codimension one surface in the energy surface (that is a 2-dimensional surface like our DS) is simply given by the interior product of the energy surface volume form *η* with the Hamiltonian vector field,

To compute *φ*, we take two vectors *ξ*_{1} and *ξ*_{2} tangent to the DS and a vector *ξ*_{3} such that *dH*(*ξ*_{3}) ≠ 0. First we note that the volume form Ω when applied to *ξ*_{1}, *ξ*_{2}, *ξ*_{3} and *X _{H}* gives

The last line follows since *ξ*_{1} and *ξ*_{2} are tangent vectors to the DS and therefore *dH*(*ξ _{i}*) = 0 for

*i*= 1, 2. Also, we have

*dH*(

*X*) = 0 as

_{H}*X*is a tangent vector field to the energy surface. From the definition of the energy surface volume form, we have

_{H}Since *φ*(*ξ*_{1}, *ξ*_{2}) = *i*_{XH}*η*(*ξ*_{1}, *ξ*_{2}) = *η*(*ξ*_{1}, *ξ*_{2}, *X _{H}*), it then follows from (B8) and (B9) that

and

The flux form through the DS is simply the symplectic 2-form applied to tangent vectors of the DS.

*Condition for non-tangency of the Hamiltonian vector field*. To check that the vector field is not tangent to the DS (i.e., trajectories cross the DS) we need to evaluate the flux form on tangent vectors to the DS.

The calculation we need to carry out is the following:

The flux through the DS is then zero whenever

That is, the ratio of velocities corresponds precisely to the phase point being on the NHIM (PO).

## REFERENCES

Consider a potential energy function *V* = *V*(*q*_{1}, …, *q _{n}*) that is a function of

*n*coordinates {

*q*}. (Coordinates describing translation and rotation are excluded.) At a non-degenerate critical point of

_{k}*V*, where ∂

*V*/∂

*q*= 0,

_{k}*k*= 1, …,

*n*, the Hessian matrix ∂

^{2}

*V*/∂

*q*∂

_{i}*q*has

_{j}*n*nonzero eigenvalues. The

*index*of the critical point is the number of negative eigenvalues.

A periodic orbit is said to be neutrally stable if all of its Lyapunov exponents are zero.