Identifying Reaction Pathways in Phase Space via Asymptotic Trajectories
Yutaka Nagahata, Florentino Borondo, Rosa M. Benito, Rigoberto Hernandez
PPhys. Chem. Chem. Phys.
Identifying Reaction Pathways in Phase Space via Asymptotic Trajectories
Yutaka Nagahata, F. Borondo,
2, 3
R. M. Benito, and Rigoberto Hernandez ∗ Department of Chemistry, Johns Hopkins University, Baltimore, MD 21218 Instituto de Ciencias Matem´aticas (ICMAT), Cantoblanco, 28049 Madrid, Spain Departamento de Qu´ımica, Universidad Aut´onoma de Madrid, Cantoblanco, 28049 Madrid, Spain Grupo de Sistemas Complejos, Escuela T´ecnica Superior de Ingenier´ıa Agron´omica,Alimentaria y de Biosistemas, Universidad Polit´ecnica de Madrid, 28040 Madrid, Spain (Dated: January 6, 2021)
Abstract
In this paper, we revisit the concepts of the reactivity map and the reactivity bands as an alternative tothe use of perturbation theory for the determination of the phase space geometry of chemical reactions. Weintroduce a reformulated metric, called the asymptotic trajectory indicator, and an efficient algorithm toobtain reactivity boundaries. We demonstrate that this method has sufficient accuracy to reproduce phasespace structures such as turnstiles for a 1D model of the isomerization of ketene in an external field. Theasymptotic trajectory indicator can be applied to higher dimensional systems coupled to Langevin baths aswe demonstrate for a 3D model of the isomerization of ketene.
Keywords: transition state theory; phase space geometry; normally hyperbolic invariant manifold (NHIM); reactivitymap; reactivity bands; Langevin Equation a r X i v : . [ phy s i c s . c h e m - ph ] J a n . INTRODUCTION The identification of a reaction path (or pathway) has received attention from the beginning ofthe development of transition state theory (TST) to characterize the energetics of the reactionbetween reactants and products. Eyring called it “the path requiring least energy” which is nowcommonly called the minimum energy path (MEP) and obtained as the path of least resistancestarting from the energy minimum associated with the reactant. Fukui is generally credited forwhat came to be known as the intrinsic reaction coordinate (IRC) because he introduced it, in amass-weighted coordinate system, as the path of steepest descent starting from the saddle. Beyondthe developments in the statistical formulation of TST, the re-imagination of the reactionpath as an object in full phase space led to the use of reactivity bands and the periodicorbit dividing surface (PODS) to characterize reactions. Both analyses are formulated in termsof identifying the trajectories in between the reactive and non-reactive trajectories. See Ref. 15for more details, and Ref. 16 for the connection to reaction path sampling also qualitativelydescribed in our earlier work. We are thus led to use these mathematical structures to reconsiderthe determination of the optimal reaction path in phase space.In the 1980s, dynamical system theory was advanced through the use of the Poincar´e map, turnstile structures, and reaction island theory. Unfortunately, all of these methods are applica-ble mostly to systems up to 2 degrees of freedom (DoF). Going beyond this restriction, Wiggins suggested a multidimensional generalization of the unstable periodic orbit in 2 DoF systems andthe reaction pathway associated with it. The former is a normally hyperbolic invariant mani-fold (NHIM), and the latter is the boundary of the reaction pathway which can be understood asthe stable and unstable manifolds of the NHIM. In the limit of two dimensional systems, the NHIMis an unstable periodic orbit which was later understood as an anchor of the PODS. Fenichel was the first to prove that these NHIMs persist under perturbation. For simplicity, in this work,we call this and its subsequent generalizations, the NHIM Persistence Theorem (see Appendix A).A consequence of this theorem is that the saddle point on the potential energy surface plays asignificant role in many cases because the structure of the NHIM near the energy of the saddlepoint will persist for larger energies as long as the truncated higher-order terms in perturbationtheory are small enough.The pioneering work on the application of perturbation theory for reaction dynamics in thechemical physics community was made by Hernandez and Miller through Van Vleck semiclas-sical perturbation theory and by Komatsuzaki et al. through classical Lie-canonical perturbation2heory (CPT) (also known as a component of Birkhoff normal form theory) to obtain non-recrossingdividing surfaces in many DoFs. Uzer et al. later showed the relation between this normal formtheory and the geometry of the reaction pathway. The normal form theory was also further general-ized to address related challenges in quantum systems, the effects under rotational coupling, Langevin dynamics, generalized Langevin dynamics, and the classical and the quantum dynamics under an external field. The extraction of the transition state (TS) trajectory andrelaxation of the normal form theory plays a crucial role especially in time-dependent and stochastic theories. Naturally, perturbation theories are limited and are applicableonly when the zeroth order approximation (typically normal mode Hamiltonian) is valid, and theasymptotic series returns converged results.Challenges to this theory arise when the TS or the phase space bottleneck is not stronglydominated by a potential energy saddle point. Perturbation theories expanded around such a pointare not expected to be accurate unless the bottleneck happens to remain within the convergenceradius of the perturbation. One possible such challenge comes from the existing of the roamingreaction pathway observed experimentally in formaldehyde H CO decomposition path to H +CO. The phase space manifestation of this system shows the significant role of the reactivityboundary in the absence of a potential energy saddle. Nevertheless, the structure of transition statetheory can be preserved even when such roaming reactions are present through the identificationand use of a global transition state dividing surfaces. There are other known examples that thereaction mechanisms are irregular due to factors such as long-time trapping in a well, bifurcationof the PODS, dynamical switching of the reaction coordinate, and the presence of near higher-index saddles, with chemical species shown in the references. In the limit of two or fewer DoFs,there are nonperturbative approaches such as periodic orbit analysis, which is used in an effectively2 DoF system of the roaming reaction path: H CO −−→ OC–––H −−→ H + CO. However, theextension of these approaches to the systems with higher than 2 DoF is still a challenging task.An alternative approach to determining the reactivity boundary is rooted in the Lagrangiancoherent structure (LCS) of the dynamics in the Lagrangian frame. The LCS is mediatedby the NHIM and its stable/unstable manifolds, and can be evaluated by numerical analysis — e.g. , through the identification of the finite time Lyapunov exponent (FTLE) ridge. The FTLEanalysis has mostly been employed in effectively two-dimensional systems, such as ocean flows dueto theoretical and numerical limitations. The Lagrangian descriptor (LD) was proposed as aheuristic alternative. In the context of reaction dynamics, the LD was first introduced by Cravenand Hernandez. The theory provides, for example, a way to obtain nonrecrossing dividing surfaces3n barrierless reactions in which perturbation theory is nonsensical because no unique zeroth-orderdividing surface is available. It has also been used to reveal geometric features in molecularsystems such as ketene and LiCN. In this paper, we propose a non-heuristic approach for locating the reaction pathways in a phasespace. The approach is an alternative to the nonperturbative methods cited above. Toward itsformulation, we revisit the theory of the reactivity map and reactivity boundaries. We introducethe asymptotic trajectory indicator (ATI) associated with the reactivity map, and formulate itin the context of dynamical systems theory in Sec. II. It is employed in a numerically efficientalgorithm to extract the reactivity boundary in Sec. III. The analysis is applicable whenever thesolutions of the equations of motion are continuous in some sense with respect to initial conditionsand integration time. In the case of smooth Hamiltonians, this is automatically satisfied withrespect to the standard definition of continuity, but in more general cases, such as with stochasticequations of motion, it is sufficient to define continuity with respect to neighborhoods of the inputand output variables. Consequently, one can use ATI even for systems coupled to a Langevin bathas discussed briefly in Subsec. II E.In what follows, we demonstrate the neighbor bisection and continuation with ATI (NBC-ATI)method through application to the 1 DoF and 3 DoFs reduced ketene models under an externalfield and a Langevin bath, respectively, in Sec. IV. In the former case illustrated in Fig. 1a, thismethod uses the ATI shown in Fig. 1b to uncover the complex reaction pathway —with respectto the reactivity boundary— and how it is guided by the phase space structure associated withthe four potential energy saddles shown in Fig. 1a. As a result of this analysis, we find the rarereactive pathway between ketene and the intermediate structures —that is, formylmethylenes andoxirene— shown as the black dashed line in Fig. 1c. In general, the ATI can be used to locate thephase space skeleton of the reaction, the stable and unstable manifolds of all available NHIMs. II. THEORYA. Phase Space Flow around an Index-one Saddle
For the normal mode approximation of a Hamiltonian expanded at an index-one saddle point: H ( p , q ) = 12 ( p − ω q ) + 12 n (cid:88) i =2 ( p i + ω i q i ) , (1)here we use i th normal mode coordinate q i , its conjugate momentum p i , its frequency ω i , thecoordinates vector q = ( q , . . . , q n ), and the momenta vector p = ( p , . . . , p n ). For all i , ω i is4 xirene CC HH O formylmethylene CC HHO formylmethylene C HH C O ketene C HH O C ketene CHH OC π /3 π /30- π -2 π /3- π /3 π min1 min2 min3 FIG. 1. Highlights of the main results of this paper for the 1 D ketene isomerization reaction in an externalfield: (a) The potential energy V F ( q F ) (solid black line), the external force through dipole interaction V ex ( q F , t ) at phases ωt = nπ/ n = 3 , , . . . , − V ex shown in the rightaxis), the TSs (gray plus symbols), the minima (gray filled-circles and vertical dashed lines) corresponding tothe chemical formulas (ketene, formylmethylene, and oxirene), and the left (red) and right (blue) absorbingboundaries. (b) The first passage time to the left (right) absorbing boundary τ Lf ( τ Rf ), or recurrent timeafter ωt = 14 π (gray), or ATI with the color scale given at the bottom. They provide the phase spaceskeleton of the reaction (yellow). (c) The phase space geometry structures shown here for comparison arethe stable (blue) and unstable (red) manifolds, and regions whose points exit right (blue) or left (red). Asample trajectory (black dashed line) is shown to illustrate the points on the Poincar´e surface of section incorresponding regions: square and circles on the time slice | ωt | ≡ π ). positive so that ω i is a pure imaginary frequency. The index of a potential energy saddle point5 b ) -1-2 0 1 2 q -2-1012 p ( a ) E : E n e r gy o f r eac ti v e D o F qp FIG. 2. The phase space flow of trajectories along a reactive DoF is shown for various energies and itsprojection in phase space. For simplicity, the Hamiltonian, H = ( p − q ) /
2, of this hyperbolic normalmode is taken to be that at second order in the potential expansion. Reactive trajectories (red, E = 1),non-reactive trajectories (blue, E = − E = 0) in (a) the q - p - E space, and (b) q - p space. The potential energy surface is also drawn in gray on (a) q - E plane. is determined by the number of Hessian eigenvalues, which correspond to the Morse index of acritical point. The Hamiltonian equation of motion can be written with ( p , q ) bydd t p i q i = ∓ ω i q i p i , (2)6here ∓ correspond to the sign of the monomial ± ω i q i in Eq. (1), respectively. The first normalmode i = 1 is called hyperbolic or reactive DoF, and is shown in Fig. 2 with ω = 1. In the upperpanel (2a), trajectories with given reactive mode-energies are shown relative to the correspondingpotential energy saddle. When the mode-energy is positive and negative, the trajectory is reactiveand non-reactive, respectively. The boundaries in between consist of asymptotic trajectories withzero mode-energy toward or from the origin. Here, the origin of the reactive DoF is called theNHIM, and its dimensionality grows as the total number of DoFs is increased. The trajectoriesasymptotic toward or from the NHIM are called the stable or unstable manifolds of the NHIM,respectively. In Wigner’s TST formulation, the dividing surface is defined as q = 0 and p > p , q ) = (0 ,
0) is known as an anchor of the dividing surface. For a normalmode Hamiltonian, one can solve the equation of motion by the finding constants of motion, i.e.,normal mode action J . One can rewrite the reference Hamiltonian as H = n (cid:88) i =1 ω i J i ( p i , q i ) , (3)where the i th generalized momentum is J i = ω i ( p i − ω i q i ) ( i = 1) ω i ( p i + ω i q i ) ( i ≥ . (4)which is conjugate to the action θ i , and together satisfydd t J i θ i = − ∂H / ∂θ i ∂H / ∂J i = ω i (5)with its solution θ i = ω i t, J i = constant . .If H is dominant, a similar relation can be obtained from CPT. When one can expand theHamiltonian at the index-one saddle point: H ( p , q ) = H ( p , q ) + (cid:88) k =1 (cid:15) k V k ( q ) , (6)where V k ( q ) is a ( k + 2)-order polynomial in q , and the perturbation order k is tracked by (cid:15) = 1without changing the equation.The construction of perturbation theory now follows a series of canonical transformations thatsuccessively remove the terms in θ up to a desired order while formally preserving the Hamiltonianstructure. That is, we seek to find the composite transformation from ( ˆ p , ˆ q ) to ( ˆ J , ˆ θ ) such that thenew Hamiltonian is H = ˇ H ( ˆ J )+ O (cid:0) (cid:15) k (cid:1) . In Lie-CPT, this is achieved through a “time” propagation7 IG. 3. Forward extremal Lagrangian Descriptor L f (arc length of a trajectory) over initial conditions of a2-DoF harmonic system: H = (cid:80) i =1 , p i / (2 m i ) + m i k i q i / m = 0 . , k = 0 . , m = 5 , k = 3, with E = 10, q (0) = 0, and p (0) >
0. The light-yellow area is energetically prohibited. (a) time evolution of L f ( t ) for q (0) = 0 in t − p (0) space, (b) L f ( t = 1000) in q (0) − p (0) space. that go forward or backward in time resulting in the solutions, ˆ F and ˇ F , respectively. One cansolve the equation of motion of ˇ H with the order of accuracy O (cid:0) (cid:15) k (cid:1) , i.e.,dd t ˆ J i ˆ θ i = − ∂ ˇ H (cid:46) ∂ ˆ θ i ∂ ˇ H (cid:46) ∂ ˆ J i = ω i , (7)where ˜ ω i := ∂ ˇ H (cid:46) ∂ ˆ J i . This equation can be rewritten with ( ˆ p , ˆ q ) asdd t ˆ p i ˆ q i = ˜ ω i ( ˆ J ) ω i ∓ ω i ˆ q i ˆ p i , (8)where the difference from Eq. (2) is the locally constant term, ˜ ω i ( ˆ J ) /ω i . Therefore, the coordinatetransformation ( p , q ) → ( ˆ p , ˆ q ) gives a local (in the sense of O (cid:0) (cid:15) k (cid:1) ) independence for each DoFincluding the reactive DoF ( i = 1). This suggests that the phase space flow of Eq. (8) has asimilar shape with Fig. 2, but in the space of ( ˆ p , ˆ q ). Revisiting the flow of the trajectories withoutspecifically constructing the Lie-CPT transformations should thus result in an alternate construtionrevealing the reaction path, and serves to motivate the approach pursued here. B. Inhomogeneity of the LD on the NHIM
To obtain the NHIM and its stable and unstable manifolds non-perturbatively, the LCS and LDare introduced here. In dynamical systems theory, a distinguished hyperbolic trajectory has beendefined as the non-autonomous analogue of a hyperbolic fixed point. Mancho and coworkers implemented the extremal LD , M ex , which describes the arc length of trajectories in coordinate space M ex ( q , ˙ q , t ; t ) := (cid:90) t + tt − t (cid:107) ˙ q (cid:107) d t . (9)The arc length also evaluated for the forward and backward LDs: L f := (cid:90) t + τt (cid:107) ˙ q (cid:107) d t (10a) L b := (cid:90) t t − τ (cid:107) ˙ q (cid:107) d t (10b)to estimate stable and unstable manifolds, respectively, as the “abrupt change” of LDs.The LDs, L f and L b , —refer to Eq. 10— are accumulated value of a positive scalar along atrajectory. Trajectories that diverge from each other will necessarily accumulate different LDs, andthe separation of these values appear to signal the presence of a stable/unstable manifold of theNHIM that lies between them. Although there are known examples that the singular contour ofthe LD correctly corresponds to the stable/unstable manifolds, there is no general proof on thecorrespondence. Formally, the LDs would be obtained at t → ∞ where the values all go to infinityregardless of the choice of trajectory. The exceptions arise from fixed points at which the LD oftrajectories asymptotic to them converges to a finite value. In practice, we integrate for a long, butfinite, time at which there is a visible feature, deviation, or “abrupt change” in the LD betweeninitially nearby trajectories, or abrupt features in the LD such as narrow ridges or valleys.There are some concerns about the generality of the conjecture. In particular, Haller criticizedthe use of the LD because it is not objective. That is, as long as the LD is defined by a norm,such as an arc length, the LD value is not necessarily independent of the particular choice of DoFs.To illustrate this concern, let us consider a three-dimensional normal mode Hamiltonian whichhas one reactive DoF and two vibrational DoFs. If the trajectory is on the NHIM, the dynamicalvariables ( q, p ) of the reactive DoF remain constant, i.e. , they remain on the saddle with zeroreactive velocity. If LD values are uniform over the two vibrational DoFs, the reactive DoF is thenthe only relevant DoF. However, this is not the case as shown in Fig. 3. On the other hand, forthe reactive mode ( q , p ) of a normal mode Hamiltonian, a trajectory on the stable and unstablemanifolds of the NHIM, with a momentum expressed by p (0)e − λt and p (0)e λt , reaches the NHIMat t → ∞ and t → −∞ respectively. This means that around the NHIM, extremal LD can beaffected by the reactive DoF and be prone to dominant contributions from the vibrational DoFs.9 IG. 4. The forward ATI, τ R , Lf (left column), and backward ATI, τ R , Lb (right column), defined in Eqs. (16)and (17) respectively. They correspond respectively to the first hitting time of the right q = 3 (blue) andleft q = − In some cases, modification of some initial conditions along the manifold can have larger effects onthe LD value than from those not along the manifold. For these cases, the ‘abrupt’ change in theLD value could misidentify the region as containing a NHIM. Thus the use of the LD to identifyNHIMs has to be done with care.
C. Asymptotic Trajectory Indicator
An alternative to the LD can be achieved through the reactivity bands (in one-dimensionaldomains) or reactivity map (in higher dimension) and reactivity boundary. These struc-tures are defined on the domain of initial condition in phase space by designating them accordingto their ultimate origin or destination to a reactant domain or one of possibly many distinct prod-uct domains. Between initial conditions assigned to different final basins, there could be initialconditions whose trajectories never reach one of the designated basins, and thus act as reactivityboundary. For example, the purple and pink trajectories in Fig. 2 form the boundary betweeninitial conditions assigned to products (at t → + ∞ ) and reactants (at t → −∞ ), respectively. To-gether, these boundaries separate trajectories into four categories whose initial and final chemicalstate can be assigned according to whether trajectories are (1) staying inside of, (2) entering into,(3) exiting from, or (4) staying outside of the trapping area. Those structures are fundamental tothe turnstile and reaction island theories. 10o compare with the NHIM theory, let us formulate the reactivity boundaries mathemat-ically by using their asymptotic nature. Let φ t be the propagator in time t , i.e. , a flow function , φ t : x ( t ) → x ( t + t ) . (11)For given φ t , there could be a set of trajectories, which is restricted to a subspace M . Such asubspace M of the phase space P is said to be invariant when M := (cid:8) x (cid:12)(cid:12) ∀ x ∈ M and ∀ t, φ t ( x ) ∈ M (cid:9) . (12)Or more simply, M is an invariant subspace when φ t ( M ) = M for arbitrary t . If there areasymptotic trajectories to M then one can define the stable manifold and unstable manifold of theinvariant manifold M as follows: W (s) M := (cid:8) x (cid:12)(cid:12) ∀ x ∈ P , lim t →∞ φ t ( x ) ∈ M (cid:9) , (13) W (u) M := (cid:8) x (cid:12)(cid:12) ∀ x ∈ P , lim t →−∞ φ t ( x ) ∈ M (cid:9) . (14)Accordingly, W (s) M ∩ W (u) M is invariant. For simplicity, we define the reactivity boundary separatingthe destination and origin of trajectories as W (s)asym and W (u)asym respectively. This allows us todesignate an invariant manifold related only to W (s)asym and W (u)asym as M asym := W (s)asym ∩ W (u)asym , (15)where W (s) X := W (s) M X , and W (u) X := W (u) M X for X = asym.In practice, this manifold M asym can be detected through observation of nearby trajectories.We first identify a subspace S that contains M asym — i.e. , S ⊃ M asym ,— with initial conditionsassociated with trajectories that remain in W (s)asym under forward propagation or W (u)asym underbackward propagation. We then define the first passage time for each point x ∈ S according towhen it first crosses the boundary ∂ S , that is, τ f ( x ; S ) := min t (cid:8) t (cid:12)(cid:12) ∀ t ≥ , φ t ( x ) ∈ ∂ S (cid:9) . (16)Consequently, for x ∈ W (s)asym , then τ f ( x ; S ) = ∞ . Neighbors y of x ∈ W (s)asym , will necessarilyhave large but finite τ f ( y ; S ) because of the continuity of the equation of motion. Therefore,abrupt changes in τ f ( x ; S ) indicate the nearby location of the singular contour , i.e. , W fsing := { x | τ f ( x ; S ) = ∞} . This contour is formed from forward asymptotic trajectories, and is thereforethe manifold W ( s )asym . For this reason, hereafter we call τ f ( x ; S ) the forward asymptotic trajectoryindicator (ATI). Similarly, we can define τ b ( x ; S ) := max t (cid:8) t (cid:12)(cid:12) ∀ t ≤ , φ t ( x ) ∈ ∂ S (cid:9) (17)11s the backward ATI . We provide an example in Fig. 4 to show a typical behavior of the ATI.It is useful to consider how the ATI generalizes for more complex cases, including those when ∂S encloses multiple asymptotic manifolds. In the following, we prove that when one observesan ATI value of a point on ∂S , its value is almost surely finite, if S is compact. Here, “almostsurely” means that the statement is true except for one or more dimensions less than an equi-energysurface of the phase space. First, let us consider a subset of initial conditions on ∂S , for which entiretrajectories are not bounded to the region S . Each such trajectory must have an even number ofintersections on ∂S —that is, they go in and out in pairs— because it is not bounded. Second,the other initial conditions are known to be almost surely recurrent as specified by Poincar´e’srecurrence theorem . (See, e.g. , Sec. 3 16 D in Ref. 69.) Asymptotic trajectories serves as examplesof such measure zero sets. Therefore, when one observes an ATI value of a point on ∂S , its value isalmost surely finite, if S is compact. In fact, below we observe the reactivity boundaries as singularcontours even if S includes multiple asymptotic manifolds. In practice, the finiteness of the ATImay indicate the requirement of very long time propagation. For example, a trajectory may betrapped in a potential energy well. The complexity of highly coupled reactions also challenges ourapproach because they can lead to fractal structures arising from chaotic dynamics. Such casescan be avoided by taking S around a M asym . In Subsec. IV B, we further show heuristics to avoidlong-time trajectory calculations. D. Comparison with NHIM theory
Here we confirm that the ATI leads to a NHIM in cases when the solution is accessible to pertur-bation theory, and how it extends beyond it. To this end, we first observe that the NHIM PersistenceTheorem provides for the persistence of the phase space geometry —vis-`a-vis the NHIM— underperturbation. As we reconfirm in the Appendix A, the reaction dynamics around an index-onesaddle on a potential energy surface satisfies the requisite conditions of the theorem. In this case,the NHIM corresponds to a DoFs orthogonal to the reactive DoF at the TS. The latter is thenonlinear analogue of the DoF associated with the eigenvector of the positive Hessian eigenvalue.The stable and unstable manifolds associated to the reaction coordinate are illustrated in Fig. 2.They are associated with imaginary frequencies, called characteristic exponents throughout thiswork, that characterize their strongest decay. That is, the persistence of the manifolds resultsfrom their exponential expansion and contraction, respectively, and the sign of the characteristicexponent reflects this. Suppose for its imaginary frequency λ s , u and upper bound of the other12maginary part of frequencies λ c , the NHIM has r ≥ ≤ rλ c < λ s , u ( e.g. λ c = 0 whenthe saddle is index one). In this case and if the flow is k -differentiable φ t ∈ C k , then M NHIM ∈ C k when k ≤ r . In addition, for a f ∈ C k , f ( M NHIM ) is still a NHIM and k -differentiable, i.e. , persistent under C k perturbation. Specifically, this NHIM is called a k -NHIM. (SeeAppendix A for the explicit definition).For example, if M ( ⊂ S ) is a NHIM associated with an index-one saddle, the characteristicexponents − λ s , u and λ s , u corresponding to expansion and contraction directions are associatedwith stable W (s)NHIM and unstable W (u)NHIM manifolds, respectively. The singular contour of the ATIsdefined in the previous section is a stable ( i.e. , W fsing ) or unstable ( i.e. , W bsing ) manifold of theNHIM.More generally, M NHIM and M asym are not equivalent for higher index saddles. For example,if there are other exponentially growing directions with characteristic exponents − µ s and µ u inaddition to those of that stable − λ s and unstable λ u manifolds, satisfying − λ s < − µ s < < µ u <λ u , then the directions are tangent to M NHIM because of λ c ≥ max { µ s , µ u } . On the other hand,those directions are still normal to M asym because they converge to it asymptotically. Therefore, M asym is not a NHIM, and M asym may not have persistence under perturbation. Such a failurewas recently seen for an index-2 saddle despite the convergence of perturbation theory. The reactivity boundary ( W ( x )asym with x = s , u) is relevant to reaction dynamics because itmarks the transition at which the fate of the reactants is cast. For example, for a reaction with anindex-1 saddle, the reactivity boundaries are the same as the stable and unstable manifolds of theNHIM. This feature allows us to detect the boundary of reactivity even for reactions associatedwith higher index saddles. E. NHIM of Random Dynamical Systems
The nature of a trajectory resulting from a stochastic differential equation (SDE) is completelydifferent from an ordinary differential equation (ODE).
For example, a particle whose motionis described by a Langevin equation is thereby driven stochastically by a random force ξ ( t ) onlyin a formal sense. In practice, one must integrate over the latter to propagate the particle withina difference equation, e.g. , by the Euler-Maruyama method. Integrals of the random force areneither smooth nor continuous in the usual sense. Its statistics generally satisfy that of a Wienerprocess. Specifically, we can define the difference in the accumulated random force across a time13nterval as B s + t − B s ≡ (cid:90) s + ts ξ ( τ )d τ ∼ N (0 , Dt ) , (18)where D is the diffusion constant. The relation X ∼ N ( µ, σ ) is defined such that X is a randomvariable resulting from a normal distribution with mean µ and variance σ . With this construction, B t has continuity in a weak sense; that is, it satisfies the β -H¨older continuity for all β satisfying β < / Recall that β -H¨older continuity requires the existence of a constant c > β > | B t − B s | ≤ c | t − s | β . (19)for all times t and s . We further require β = 1 so that the function is differentiable. In the presentcase, the existence of weak continuity allows us to claim uniqueness of the solution for a givenstochastic sequence, and hence each solution has pathwise uniqueness. Second, as a consequenceof the lack of smoothness in B t , it is nowhere differentiable along t , and it can not be inverted toa ξ .The NHIM Persistence Theorem describing the geometry of the solutions of the SDE can beframed through an analysis of the paths (trajectories) of the SDE. For each d -dimensional accumu-lated random force B t , one can uniquely obtain a d -dimensional, t -continuous function ω ( t ) that isassociated with the saddle point (precisely defined in Appendix C) at each instance of time t , andfor which we are free to initialize at ω (0) = 0. (Note that this ω ( t ) is simply an abstraction of theso-called TS trajectory. ) Then a probability is determined from a bundle of instances ω . Let usalso introduce a t -origin shift to the probability measure, the so-called Wiener shift θ s , such that,( θ s ω )( t ) = ω ( s + t ) − ω ( s ) , s > . (20)This θ s introduces a shift in the initial time of a stochastic process to s . For example, any givenBrownian motion can be written as a particular manifestation of ω , such that, for B t = B [ ω ( t )], B [( θ s ω )( t )] = B [ ω ( s + t ) − ω ( s )]= B [ ω ( t + s )] − B [ ω ( s )] . (21)The relation can be found in Eq. (6.46) of Ref. 72. (See also Eq. (18) for a SDE.) We can thendefine stochastic analogues of the flow function and invariant manifold associated with an ODE tothe analogues of the a SDE: the stochastic cocycle φ t ω , φ tω : x ( s ) → x ( t + s ; θ s ω , x ( s )) , (22)14nd the random invariant manifold M ( ω ), φ t ω ( M ( ω )) = M ( θ t ω ) , (23)respectively. In Eq. (22), the stochastic path x ( · ) starts at s requiring a shift in the time originof ω ( t ) to s , and hence the need for the term θ s ω in the argument of x . A random dynamicalsystem (RDS) defined by Eq. (22) is uniquely obtained from a random differential equation (RDE) with the vector field f ˙ x = f ( θ t ω , x ) , (24)as usually obtained for ODE. One can obtain this RDS for some classes of SDEs, including Langevintype equations. (See Appendix C.)The NHIM Persistence Theorem for a RDS holds for random invariant manifolds definedby Eq. (23). The remarkable result of the theorem is that the NHIM is persistent under C perturbations and is C k -smooth at time t when φ t ω ∈ C k . The smoothness of φ t ω does not indicatethat the variables of SDE — e.g. of the Langevin equation— are smooth over integration time. Onthe other hand, SDE still has a H¨older continuity as detailed above. Because of this continuity,the closer an initial condition is to the stable or unstable manifold of a NHIM, then the longer thetrajectory will spend in time around the NHIM. In this way, there can still exist abrupt changesand singular contours in the ATIs of a SDE, where the latter corresponds to the NHIM and itsstable and unstable manifolds. Because of the theorem, these manifolds are smooth.As was suggested in Refs. 68 and 74, the NHIM Persistence Theorem holds for autonomoussystems. For such systems, there is a single ω which corresponds to the time series of the externalforce and does not require any fundamental change to satisfy the conditions needed to satisfythe NHIM Persistence Theorem. III. NEIGHBOR BISECTION AND CONTINUATION WITH ATI
Thus far, we have considered the case of a single invariant manifold M asym defined in Eq. (15)that lives within the subspace S of the phase space P without considering the computationalrequirements for its implementation. In practice, the latter is exacerbated by the existence ofmultiple NHIMs in S . Nevertheless, typically some features of its structure (such as possiblelocation or locations) are approximately known, and they can be used to optimize sampling ofcandidate points of the manifold, and improve the numerical implementation of the search. Here,15 b ) forward ( c ) backward t t T T t ( a ) FIG. 5. A schema of possible trajectories from an initial point in the region relative to the absorbingboundary surfaces is shown in panel (a) as a function of time: A trajectory hitting the left (in red) or right(in blue) absorbing boundary, and a recurrent trajectory (in gray) which returns to the same coordinatewith the same sign in the velocity. The projection of these trajectories in forward and backward time ontothe phase space is shown in panels (b) and (c), respectively. we present a practical approach for defining the absorbing boundaries, sampling the surroundingneighborhood of the reactivity boundary, and visualizing the ATI. We introduce the NBC-ATImethod to effectively increase the resolution and accuracy of the reactivity boundaries.
A. Absorption and Visualization
The values of the forward τ f and backward τ b times defined in Sec. II depend on the initialand absorbing conditions. However, the location of the singular contour in τ corresponding tothe reactivity boundary is independent of the absorbing conditions, as we described in Sec. II C. Anaive absorbing condition can be defined through a coordinate q at a value that is significant to thedynamics ( e.g. , at a potential energy minimum) and sufficiently far away from the TS. Figure 5illustrates an example in which trajectories are absorbed and thereafter assigned a first hittingtime τ —viz, the ATI.The domain of initial conditions can be labeled using a color that denotes one of the surfacesshown in Fig. 5 to which they absorb, and an intensity commensurate with the value of the ATI, τ .The result for a normal mode Hamiltonian is illustrated in Fig. 4. Initial conditions on the domainare labeled in red or blue according to whether trajectories are absorbed on the left at q = − q = 3, respectively. The intensity of the color is commensurate with the value of the ATI, τ Lf or τ Rf . In the context of the theory described in Sec. II, the absorbing boundaries at R and Lemployed here are examples of the two disconnected subsets of ∂ S .16 f =2 (b) (c) (d)(e) (f) (g) n b =2 n b =1 n b =3 n f =1 n f =3 (a) (h) n b =18 FIG. 6. Schema of selected steps in the sampling algorithm as described in the text. In any given step,the active straddling pairs (red), points (red) and area (magenta), are highlighted in color, and previouslysampled straddling pairs (gray), area (light gray) and point (gray) are highlighted in gray scale.
As shown in Fig. 4a, τ f or τ b has a singular contour for the forward or backward asymptotictrajectories that corresponds to the purple or pink trajectories in Fig. 2, respectively. Hence, theabrupt change in τ f or τ b in Fig. 4b indicates that there is, indeed, a reactivity boundary nearby.This singular contour is the stable or unstable manifold of the NHIM. For the present case, theNHIM is ( q, p ) = (0 , p = − q , and the unstable manifold is p = + q . B. Efficient Sampling Algorithm
When analytical methods, such as perturbation theory, fail to produce the exact form of thereactivity boundary, we must resort to using numerical methods. The challenge to the numerics,however, lies in the fact that the reactivity boundary has measure zero, and hence statisticalsampling is inefficient. To overcome this challenge, we employ an algorithm which effectivelygeneralizes the one-dimensional bisection method to 2D and the required higher dimensionalityof the space that contains the reactivity boundary. It may be implemented iteratively when auser needs to improve resolution to, for example, increase the number of points of the reactivityboundary.To initiate the first stage of the algorithm, we need a low-resolution representation of thereactivity boundary. We can construct this by way of creating a low-resolution grid in the domain,and performing a “brute-force search” across all the vertices to identify pairs of points that straddlethe reactivity boundary in the sense that one of the two adjacent points goes to the reactant-sidesurface and the other to the product-side surface. We call such pairs of points straddling pair s.In more detail, we execute the brute-force search on a (2 N BF + 1) × (2 N BF + 1) grid ( N BF = 1in Fig. 6a). The grid results in a spacing resolution at L / N BF where for the length of given 2Drectangle domain L = ( L x , L y ). As long as the reactivity boundary does not fold within a width17f this resolution, a straddling pair will capture one of its points within its connecting segment.In the first stage — “seed refinement” —, we refine the initial straddling pairs up to thedesired resolution. Specifically, we iterate the bisection method N res − N BF times to find a newset of straddling pairs to the desired grid resolution, i.e. , (2 N res + 1) × (2 N res + 1). In this stage,we assume that the desired structures are all larger than the resolution of the initial grid. Theresulting resolution of the straddling pairs, illustrated in Fig. 6b-d for N res = 3, is on the order of L / N res .In the second stage — “path-following” or “continuation” —, we construct edge-to-edge andrecurrent chains of straddling pairs. To construct edge-to-edge chains, we first choose one ofstraddling pairs which are on the edge of the domain. If it exists, then we use the pair to identify anext pair (red line in Fig. 6e-h) that also straddles the reactivity boundary. A chain of straddlingpairs is constructed by repeating this iteration until it reaches the edge of the domain (Fig. 6h).If after this construction, there remains straddling points on the edge, we pick one of them andapply the path-following stage above until they are exhausted. To construct recurrent chains, wethen repeat the procedure on remaining straddling pairs inside the domain, which will necessarilyend on themselves, until all straddling pairs are exhausted.In the third stage — “precision refinement” —, we refine the precision of the straddling pairs ofthe obtained chains. We apply the bisection method to each identified straddling pair N pre − N res times to achieve the desired grid precision (2 N pre + 1) × (2 N pre + 1) with L / N pre . In the numericalapplications, we chose N pre = 30. Because 2 − ≈ − and the highest resolution in double-precision is 16 decimal digits, then the points between straddling pairs can be differentiated onlyup to 7 additional digits. If we use 10 − for the desired accuracy in the time integration, thenwe still have 4 reliable such digits in a single iteration of time propagation. Thus, for this setup,the computational result is reliable as long as the accumulated numerical error through the timepropagation for obtaining τ is less than a factor of 10 times the error of a single step iteration.Suppose that we use a (2 N BF + 1) × (2 N BF + 1) grid in defining our initial search space, and wewanted to get a resolution of δ L / N res . Naively, this would require the determination of (2 N res + 1) points on a two-dimensional grid. Using our algorithm, instead, we expect that the number of pointsis approximately proportional to 2 N res . This estimate is based on an assumption that the numberof straddling pairs on a (2 N res + 1) × (2 N res + 1) grid C st [2 N res + 1; M ] is proportional to 2 N res ,where M is the manifold whose straddling pairs we are identifying. The assumption is triviallycorrect when M is one-dimensional on the observed two-dimensional domain. In that case, theorder estimation can be made with the relation: C st [ x ; M ] ≤ C sq [ x ; M ] ≤ (9 / C st [ x ; M ], where18 ABLE I. Number of evaluated trajectories in each stage of the sampling algorithm. C st [ resolution ; M ]and C sq [ resolution ; M ] are the number from straddling pairs and that from covering squares of M in thegiven resolution , respectively.Stage Method Number of evaluated trajectoriesExact Order (expected)0th Brute-Force (2 N BF + 1) N BF C st [2 N BF + 1; M ] × ( N res − N BF ) 2 N BF × ( N res − N BF )2nd Path-Following C sq [2 N res + 1; M ] − C st [2 N BF + 1; M ] 2 N res C st [2 N res + 1; M ] × ( N pre − N res ) 2 N res × ( N pre − N res ) C sq is the number of points attach to the squares with resolution x covering M , and the first andthe second equality are from the cases when all the straddling pairs have a different and the sameorientation from that of the adjacent pair, respectively. We know that the N -fold application ofthe bisection method to a straddling pair needs the calculation on N additional points, and thenumber of input seeds for the path-following stage is C st [2 N BF + 1; M ]. We listed the order ofmagnitude for the number of points for each stage in Table I. Based on this table, if we chose2 N BF (cid:28) N res , then the order can be estimated as 2 N res . When we include the third stage in theestimation, the order becomes N pre − N res + 1 times larger than the order up to the second step.For n DoF systems at a constrained energy, the number of 2D slices required to visualize theentire reactivity boundary in phase space becomes N (2 n − − > n ≥ N is the number ofpoints sampled along a given axis. Thus, a naive estimation of the computational cost is on the orderof (2 N res ) (2 n − − if the cost of each slice is same as above: 2 N res and N = 2 N res . Not coincidentally,this exponent (=(2 n − −
1) is the same as the dimension of the reactivity boundaries, i.e., W (s)asym and W (u)asym . In Fig. D.3 of Appendix D, we illustrate how the order estimates surmised herecorrespond to the computational cost seen in practice.The output of the algorithm — i.e. , the chain of straddling pairs— can be used to sample otherparts of the boundary by propagating forward or backward in time (See Fig. D.4 for example).However, due to the nature of a trajectory on a stable or unstable manifold, the distance betweenadjacent straddling pairs will exponentially increase by the propagation. Thus, depending on theintegration time, we may need to obtain more pairs in between some pairs we already have. In suchcases, we can use the output as an input to the first stage of the algorithm. We can then applythe bisection method N (cid:48) res times father, and thereby obtain 2 N (cid:48) res times greater resolution than theinput as the output. The interactive application of our algorithm thus results in an approximate19 ABLE II. Parameters of the reduced ketene model of Eq. (25) taken from Ref. 55.Parameter Value Unit a − . × − E h a − a . × − E h a − a − . × − E h a − c . × − E h a − d . a − m F . m e m H . m e k . × − E h a − d − . × − E h a − k . × − E h a − d − . × − E h a − curve representing the reactive boundary to a desired resolution (by the second stage) and precision(by the third stage.) IV. NBC-ATI ANALYSIS OF KETENE ISOMERIZATIONA. Reduced Ketene Model
To illustrate the NBC-ATI method, we apply it to the reaction dynamics of ketene to uncover itsreaction geometry by ATIs. To this end, we use a reduced model of ketene introduced by Gezelterand Miller and adopted by Craven and Hernandez to illustrate LDs. V ( q F , q , q ) := V F ( q F ) + (cid:88) i =1 , k i (cid:18) q i + d i q k i (cid:19) , (25a) V F ( q F ) := a q + a q + a q + cq e − dq . (25b)We use the parameters fitted by Gezelter and Muller reproduced in Table II with the correctionon m F by Ulusoy et al. The fit is based on the result of an ab initio calculation at the levelof CCSD(T)/6–311G( df,p ) by Scott et al., where q F , q and q are coordinates of the systems.This model uses the normal mode coordinate of the oxirene geometry. That is, q F is the motionalong the normal mode reaction coordinate, q is the out-of-plane wagging and twisting motionof the hydrogens relative to CCO plane, and q is the in-plane rocking and scissoring motion of20 ABLE III. Parameters of the external force of Eq. (26) taken from Ref. 53Parameter Value Unit µ . ea µ ketene . ea q . a α . a − E . E h e − a − ω . E h (cid:126) − the hydrogens relative to CCO plane. All the parameters are fixed to reproduce the structures ofoxirene and formylmethylene intermediates.Ketene is known for its remarkable photo-isomerization and photo-dissociation behavior asfirst observed by Moore et al. Its unusual reaction dynamics was discussed in the context ofroaming reactions without and with the application of external forces. In the latter case,the interaction drives a dipole along q F whose moment was obtained using B3LYP/6–311+G** byCraven and Hernandez. The resulting potential interaction, and dipole moment can be writtenas V ex ( q F , t ) = E sin( ωt ) µ m ( q F ) (26a) µ m ( q F ) := µ (e − α ( q F − q ) + e − α ( q F + q ) )+ µ ketene (26b)with the parameters reproduced in Table III.Below, we first demonstrate the NBC-ATI method for the 1DoF and 3DoF ketene models underan external force. The Hamiltonian and equations of motion (EoM) for the 1DoF system is: H F = p / (2 m F ) + V F ( q F ) + V ex ( q F , t ) , (27)dd t q F p F = ∂ p F H F − ∂ q F H F , (28)where p F is the conjugate momentum of q F , and m F is the particle mass at q F . For the 3DoF21ystem, the Hamiltonian and EoM with a Langevin bath is: H = p m F + (cid:88) i =1 , p i m H + V ( q F , q , q ) , (29) d q i d p i = ∂ p i H d t − ∂ q i H d t − γ ( p i /m i )d t + d B t ( t ) , (30)where i = 1 , , F, p i the conjugate momentum, m i the mass of q i , m i = m H ( i = 1 , γ =0 . E h (cid:126) − is the friction, and d B t ( t ) ∼ N (0 , k B T γ d t ) follows a Wiener process with the Boltz-mann constant k B and the temperature T = 300 K. Here d B t ( t ) ∼ N (0 , k B T γ d t ) means that therandom variable B ( t + d t ) − B ( t ) follows normal distribution N (0 , k B T γ d t ) with average 0 andvariance k B T γ d t . By definition, d B t ( t ) satisfies the stochastic-process version of the fluctuation-dissipation theorem. Equation (30) is deterministic for a stochastic instance. That is, Eq. (30)(or generally in Itˆo process) has pathwise uniqueness of the solution. For this reason, we use thesame noise instance { d B t ( t ) } t for all the stochastic trajectories.Equation (28) is integrated numerically by the Dormand-Prince (at 5th order) method with stepsize control (for absolute and relative error is 10 − ) implemented in the C++ boost::numeric::odeint library. Numerical integration of Eq. (30) is performed by the Euler-Maruyama (at 1st order)method with d t = 0 . B. Asymptotic Trajectory Indicator
We now demonstrate the usefulness of the ATI for the 1 DoF (Eq. (28)) and 3 DoF (Eq. (30))ketene models. To this end, we use the visualization scheme explained in relation to Figs. 2 and 5.In Fig. 7, we show the result for the 1 DoF model marking each location with the value of theATI for trajectories starting at that point and ending when they reach | q F | = 3. As can be seen,for the forward time propagation, trajectories starting from initial conditions at q F = − , p F < q F = 3 , p F >
0) have the lightest red (blue) color, in Fig. 7a, because these trajectories are absorbedimmediately. Similarly, in Fig. 7b, trajectories starting at q F = − , p F > q F = 3 , p F <
0) have thelightest red (blue) color. The areas with large | p F | , not plotted in the figure, correspond to ballistictrajectories that will not be trapped in the wells. The areas which have darker reds and bluescorrespond to trajectories that bounce back at the right and left external saddles q F = 2 . − . p F = 0 due to the application of the recurrent condition, and is an artifact of the22 IG. 7. Phase space representation of the ATI values for the 1D ketene model in an external field for the(a) forward τ R , L , Cf and (b) backward τ R , L , Cb time propagation. The color for each initial condition is chosenaccording to the absorbed position: q F = − q F = 3 (blue), and the initial coordinate with the samevelocity sign (gray), as also shown in Fig.5. way we define absorption. The gray colored area results from the limitation in the propagationtime used in our calculation to be insufficiently long to resolve these areas in terms of red and blueabsorbing boundaries. For example, in the case when we apply the recurrent condition for the latertime starting at ωt = 14 π , some gray areas in Fig. 7a are colored as it is seen in Fig. 1b. We showbelow that the manifold-like structures in this area are in fact a manifold.As discussed in Sec. II, the use of the ATI to locate the stable or unstable manifold of a NHIMin a given region becomes unclear when the region contains more than one NHIM because of thechallenge in assigning asymptotic behavior to a particular origin. To address this challenge, weconsider the effect of inserting additional absorbing boundaries. For example, the results corre-sponding to the insertion of absorbing boundaries at the potential energy minima (min1, min2,min3 defined in the caption to Fig. 1 are shown in Fig. 8. The identification of the manifold resultsfrom propagation of trajectories that are 5 times smaller than the previous figure (as indicated by23 IG. 8. Phase space representation of the ATI values for the 1D ketene model in an external field for the(a) forward τ R , Lf and (b) backward τ R , Lb time propagation. The color for each initial condition is chosenaccording to the absorbed position: left (red), and right (blue). The absorbing boundaries are located at q F = −
3, the local minima (min1, min2, min3), and q F = 3 are shown in gray corresponding to the red orblue boundaries of Fig. 5 within a given region. the smaller ATI values). The increased efficiency results from the fact that we do not need to usethe recurrent condition due to the presence of the additional absorbing conditions. We are thusable to observe the manifolds that mediate internal mechanism in addition to the escaping process.Below in Subsec. IV C and Appendix D, we demonstrate the resolution of the entire phase spaceusing this approach.The ATI can be used in systems that are stochastic/non-autonomous, with multiple DoFs,and without need for an a priori reaction coordinate. To demonstrate this fact, we present avisualization of ATI for the 3 DoFs ketene model coupled to a Langevin bath. In Fig. 9a andFig. 10a, the q = q slice contour plot of the potential energy surface is shown. The white filled-circles (plus symbols) are local minima (saddles) on the slice, and the orange filled-squares (timessymbols) are projected local minima (saddles). In Fig. 9, the initial conditions are prepared on theyellow dashed line with positive out-of-slice velocity and constant energy E = 0 . . u . ( q = q = 024 IG. 9. (a) Potential energy surface for the reduced 3D ketene model on the q = q plane. The whitefilled-circle symbols (plus symbols) indicate the position of minima (saddles) on the plane, and the orangefilled-square symbols (times symbols) correspond to projected minima (saddles). The yellow dashed line isthe surface of initial conditions. The middle and the bottom panels are the forward (b), and backward (c)ATIs on the q = q = 0 plane. The absorbing boundaries are the same as in Fig. 8 except for the outermost boundaries that are now located at q F = ±
4. The light yellow areas in (b) and (c) are on the outsideof the initial energy E = 0 . . u . . and p = p > q F = ± q F = ± q F = ±
3, possibly due to the overlap of the boundarywith the NHIMs. Since we use the fixed initial energy, there is an upper and lower limit for the25elocity not seen in Fig. 8. In comparison with Fig. 8, there is no change in the timescale of theATIs in Fig. 9 because the increase of the DoFs does not affect the timescale of the motion.Although the reaction coordinate on the potential energy function is curved, we are still ableto locate the reactivity boundaries on each cell According to the result of the ATI, the initialconditions are best given by p > p <
0) for positive (negative) time integration. This fact canbe observed from the backward time propagation case, shown in the left cell of Fig. 9c. Therein, wemostly observe bounce back trajectories from the right (min1) shown by blue colors. This occurssince we take p >
0. As a consequence, trajectories have opposite velocity along p in comparisonwith the trajectories sliding down from the left external saddle. The right cell of Fig. 9c has similarbehavior due to the symmetry of the potential energy surface.The black areas in the middle right cells of Figs. 9b and 9c correspond to trajectories thathave not finished in the computation time t = 10 . Those trajectories are, in fact, those whichstay vibrating on the initial coordinate q F plane, possibly in the attraction basin of the fixedpoints/invariant manifold. There is a line in Fig. 9c (center cells) of small discontinuities which isan artifact of the absorbing boundaries. The reactivity boundary must have a singular value of τ ,otherwise it appears due to an inappropriate absorbing boundary.Beside the absorbing boundaries considered in Fig. 9, for the multiple DoFs systems, one canuse a variety of absorbing boundaries, such as a boundary transverse to the reaction path or the oneused in the reaction island theory. In Fig. 10, we present the result for the absorbing boundariesA (red) and B (blue) that enclose the left external saddle point, and are given by q = q F + 8and q = q F + 3 respectively. The initial conditions are prepared on the boundary B with q = q by using the mass-weighted coordinates along B (˜ q (cid:107) ) and orthogonal to B (˜ q ⊥ ). We define thezero axis (˜ q = 0) by q = − q F − q ⊥ . Thus,the origin ˜ q = 0 of Fig. 10 is at ( q F , q , q ) = ( − , , q (cid:107) and ˜ q ⊥ underthe condition p = p . The result for the ATI is shown in Figs. 10b and 10c. In these panels,the reactivity boundary between blue and red area can be located even in the challenging casepresented by a Langevin bath. Initial conditions with ˜ q < p < p >
0) correspondingmostly to reacting trajectories (red) in the forward or backward time propagation as seen in Fig. 10b(Fig. 10c). This is because these initial conditions are closer to the left external saddle and havevelocities ahead to or from this saddle, respectively.26
IG. 10. (a) The potential energy surface of reduced 3D ketene model on q = q plane (see Fig. 9) andsuperimposed absorbing boundaries A (red line) and B (blue line). The green dashed line corresponds to˜ q = 0. (b) forward and (c) backward ATIs of reduced 3D ketene model in a Langevin bath with initialconditions on surface B with q = q . The light yellow areas in (b) and (c) are the outside of the initialenergy E = 0 . . u . .TABLE IV. Four types of phase space regions separated by stable and unstable manifoldsstable manifoldinside outsideunstable inside (1) staying inside (2) entering intomanifold outside (3) exiting from (4) staying outside C. Turnstile and Reaction Path
We now demonstrate a way to use the manifolds obtained by NBC-ATI in Sec. III B in thecontext of turnstile or lobe dynamics. As discussed in Fig. 2, the stable and unstable manifoldsis the destination and origin dividing boundaries respectively, of the dynamics. The areas dividedby the reactivity boundaries are categorized into four types (see Table IV): trajectories that are(1) staying inside, (2) entering into, (3) exiting from, and (4) staying outside the trapping area (achemical state). Among the four types, an area which encloses the reaction pathway is categorized27
IG. 11. The turnstile (reaction pathway) of each NHIM (saddle), mediated by the manifolds obtained atphase ωt ≡ π ) of the 1D ketene model driven by an external force. The stable (blue, cyan) andunstable (red, orange) manifolds are drawn up to | ωt | = 8 π . The manifolds associated with the saddles forleft external (a), right external (b), left internal (g) and right internal (h) are shown in (c), (d), (e), (f)respectively. The colors of the manifolds become lighter, as | t | increases. Left and right panels correspondto exiting and entering trajectories, respectively, Coherent sets of these trajectories are marked in gray, withincreasingly lighter shades indicating later phases, that are enclosed by the stable and unstable manifold. into (2) or (3) and corresponds to a time slice of the pathway.In Fig. 11, the manifolds produced from each cell in Fig. 8 are shown through the superpositionof the manifolds at the phases | ωt | = 2 πM ( M = 0 , . . . ,
4) (see Appendix D). These manifoldscorrespond to the stable (blue, cyan) and the unstable (red, orange) manifolds of the NHIM ineach cell. To expose the reaction mediated by these manifolds, we draw a set of trajectoriesprepared in an enclosed area (gray). These trajectories move coherently from one enclosed areato another and are shown (stroboscopically) at ∆ t = 2 π/ω time intervals with changing strengthof its color. Darker color indicate points captured in earlier time periods. As the color becomeslighter, the set of trajectories goes out of (or into) the trapped area in the left (or right) columnin Fig. 11. Although, we illustrate trajectories in a few selected area, this should suffice to observethat the dynamics are mediated by the manifolds as expected.The manifolds in the left and right column of Fig. 11 are the same but the initial gray-coloredareas are different as they correspond to the time slice of the exiting and entering reaction paths,respectively. They illustrate the last or first few steps of the reaction pathways. One can combine28 IG. 12. The intersections of the internal and external stable/unstable manifolds of ketene 1 D withexternal field. (a) Superposition of manifolds at phase ωt ≡ π ) up to | ωt | = 14 π which partiallyshown in Fig. 11. (b) and (c) Magnifications of areas marked in (a) where we sample the initial conditionsof trajectories (plus symbols). (d) The sampled entering (pink dotted line) and exiting (green dotted line)trajectories. The plus symbols, filled-squares, and filled-circles are time-slices at the phase ωt ≡ π )of trajectories (dotted lines). these pictures to see a reaction pathway exiting into (entering from) the internal well from (to) theoutside of the observed area. Such a pathway must be explained by the intersection of the reactionpathways mediated by internal ((b) and (c)) and external ((a) and (d)) manifolds in Fig. 11. Tosee the intersection, we draw the manifolds up to the phase M = 7 in Fig. 12. We uncover twoinitial conditions that are enclosed by the both internal and external manifolds in Fig. 12b and crespectively. The resulting pink (green) trajectory shows that pathway comes into (goes out) theinternal well directly from (toward) outside. These correspond to a case in which there is a fastcooling-down (fast excitation) process due to the presence of an external field. However, from thesize of the enclosing area of the initial conditions in Fig. 12b and c, one can also observe that theamount of such coherent initial conditions are small indicating that the event is pretty rare. Thus,there is a time scale separation for transitions between internal trapping and external trappingtrajectories since the former has less energy than the latter.In Fig. 13, we interpret the trajectories of Fig 12 in the context of turnstiles. In Fig. 13a,the limit point —shown as a green filled-square— is in the area colored by the darkest gray. Astime progresses, the green filled-circles move into ever lighter colored areas in the Poincar´e map.In the first three periods, the circles from the green filled-square all moved into areas shadedwith increasingly lighter gray. Starting with the second period, the circles also move into blueareas. These later circles move into areas with increasingly lighter blue shades until they move29 IG. 13. (a) The exiting trajectories of 1D ketene shown in Fig. 12(d) superimposed on the associatedcoherent set of trajectories shown in Fig. 11 (c) (gray) and (g) (blue). (b) the same for the enteringtrajectory superimposed on the associated coherent set of trajectories shown in Fig. 11 (b) (red) and (f)(red). The color of the set of the trajectories get lighter when the time period is later. outside of the figure. Similarly, the pink trajectory starting from a different point in the phasespace experiences a series of gray and red areas in Fig. 13b upon application of the Poincar´emap. Therefore, the trajectories sampled at the points in Fig. 12 are successfully identified by theintersection of the reaction pathways mediated by internal ((b) and (c)) and external ((a) and (d))manifolds in Fig. 11.Finally, let us revisit Fig. 1. The ATI shown in Fig. 1b was computed without applying theabsorption of recurrent trajectories until | ωt | > × π . The choice for this max time is in agreementwith the number (= 7) of periodic lobes (red or blue colored areas) of the coherent trajectoriesthat we followed in Fig. 13. The coherent structures, mediated by the reactivity boundaries andrevealed by the NBC-ATI method, are clearly visible in Fig. 1b because we compute ATI valuesfor longer times than those shown Fig. 7. The yellow lines are the superpositions of the externalstable manifolds that are also shown in Fig. 12a as blue lines. These coherent structures agreewith the those shown in Fig. 1c obtained directly through the use of global manifolds. That is, thephase space skeletons are correctly extracted by the NBC-ATI algorithm. It leads to consistent andcorrect implications on the dynamics. Therefore, all the information about reactivity is extractedonly from the manifolds obtained by the NBC-ATI algorithm.30 . DISCUSSION Here, we analyze the possible use of the NBC-ATI method to describe more general chemicalreactions based on our findings from the analysis of the 1D and 3D ketene models under variouscoupling conditions. In Fig. 9 and Fig. 10, we showed a 2D slice of the ATI in the 3 DoFs phasespace. To obtain them in the full phase space, naively one needs to achieve it for all the remainingslices. However, as we showed in Sec. IV C, one can use a periodic identity of the dynamics toobtain samples by integrating across several period(s) of time. This sampling is more efficientfor autonomous systems because the phase space is identical for all time. Although this type ofidentity reduces computational costs, the minimum computation costs must be proportional to thedimension of the manifold in any numerical analysis. This is because even if we just uniformlysample a known n -dimensional manifold, a number of points proportional to the order of the power n is needed. This is a fundamental limitation of numerical sampling techniques that perturbationtheories do not suffer.Unlike in autonomous or periodic dynamical systems, there exists no natural Poincar´e mapin systems driven by aperiodic or stochastic differential equations. Hence, to see the reactivityboundaries at another time, one does not have any resource beyond the time-propagated manifold.In addition, for stochastic differential equations, the computational cost is larger because theintegrators available for such systems are not as efficient as those for an ODE. To achieve a givenresolution in the time-propagated manifold, the resolutions used within steps of the NBC-ATI mustbe chosen carefully ensuring that the computation is efficient. Due to the nature of a trajectoryon a stable or unstable manifold, the distance between adjacent straddling pairs will exponentiallyincrease by the time backward or forward propagation, respectively. Nevertheless, the weightedsampling is known to improve the efficiency. Since the NHIM and their stable and unstablemanifolds are smooth, interpolation between the straddling pairs will reduce computational costto some extent. For engineering purposes, this types of solutions can improve efficiency, although,what types of interpolation are allowed to use is still in question. Although the NHIM and its stableand unstable manifolds are of great importance, there is no guarantee that that reactivity is alwaysmediated by them. It appears in the study of a reaction associated with higher index saddles that the manifold, which is not orthogonal to the most repulsive nor attractive direction, canproduce reactivity boundaries. However, one should be careful to correctly ascribe the physicalinterpretation of these reactivity boundaries. They may not persist under perturbation. That is,a small perturbation e.g. , a small difference of potential energy surface, may introduce a large31ynamical difference. The NBC-ATI also provides a useful reference for determining which termsin perturbation theories should be retained. For example, normal form theories are known to beasymptotic series which necessarily diverge if one includes all terms in the expansion. VI. CONCLUSION
In this paper, we have presented a formulation for the ATI and an identification of reactivityboundaries based on dynamical systems theory by revisiting the reactivity map. Tothis end, we developed the NBC-ATI method which effectively requires computational resourcesthat are proportional to the dimensionality of the manifolds. We demonstrate the feasibility andefficiency of this approach on a reduced-dimensional ketene model in 1D with external field,and in 3D coupled to a Langevin bath. The NBC-ATI method can address irregular reactionswhich are not accessible to conventional perturbation theories or other types of numerical analysis.Examples include the existence of roaming pathways, bifurcation of the periodic orbit dividingsurface, dynamical switching of the reaction coordinate, and a reaction associated with higherindex saddle(s). In the reduced 1D ketene model, we obtained complex reaction paths based on turnstiles. The turnstiles consist of stable and unstable manifolds associated with the NHIMs around the fourpotential energy saddle points. We find the internal trapping area in a chaotic sea correspondingto the two formylmethylenes and oxirene conformations. Using these phase space structures, weidentified and characterized rare trajectories which are trapped by and escape from the internalwells in a short time. We thus demonstrated that our method has sufficient accuracy to reproducethe conventional analysis when it is accessible, and generalizes to reactions with more complexreaction geometry as listed above.We have shown the applicability of the NBC-ATI method for stochastic and higher dimensionalsystems through an application to a reduced 3D ketene model. We found that if one can identifyan area which includes the NHIM or an asymptotic manifold M asym as a seed of the reactivityboundaries, then the ATI still allows us to identify the reactivity boundaries in stochastic andhigher dimensional systems. Besides the advantage of reduced computational costs, the ATI allowsus to identify structures associated with the reaction path in cases which are beyond reach to theconventional perturbation theories and numerical analysis.32 NHIM W NHIM (u) W NHIM (s) E m s E m u E m c m FIG. A.1. Invariant splitting at m ∈ M NHIM : E c m ⊗ E s m ⊗ E u m and its relation with M NHIM , W (s)NHIM , and W (u)NHIM . ACKNOWLEDGMENTS
This work was partially supported by the National Science Foundation (NSF) through GrantNo. CHE-1700749. This collaboration has also benefited from support by the European Union’sHorizon 2020 Research and Innovation Program under the Marie Sklodowska-Curie Grant Agree-ment No. 734557.
Appendix A: NHIM persistence theorem
The NHIM is a multidimensional generalization of a hyperbolic fixed point such as that associ-ated with the saddle point — viz the naive TS— in a chemical reaction. Here, we recapitulate thestatement of the theorem for the persistence of the NHIM under perturbations as originally provenby Fenichel, and generalized by others. For simplicity, we call this the NHIM persistencetheorem.Suppose that in a given system, we have identified a smooth Riemannian manifold Q , a flowon Q : φ t ∈ C k ( r ≥ Q : M NHIM . This manifold, M NHIM , is aNHIM when:1. M NHIM is invariant, i.e. , φ t ( M NHIM ) = M NHIM ,2. There exists continuous splitting for ∀ m ∈ M NHIM , T m Q = E c m ⊗ E s m ⊗ E u m , (A1)33f the tangent bundle of Q at m with globally bounded, continuous associated projections: π c m , π s m , and π u m such that splitting is invariant under the linearized flow D m φ t ( E x m ) = E xD m φ t ( m ) (A2)for ∀ x ∈ E x ( x = c , s , u) where E c m = T m M NHIM and D m φ t is differential of φ t at m , e.g. the flow determined by the normal mode for the Hamiltonian systems.3. There exists constants − λ m , s < − λ m , c ≤ ≤ λ m , c < λ m , u , c m , c , c m , s , c m , u ≥
1, such thatfor the matrix D m φ t = ( ∂ x i φ tj ( x ) | x = m ) ij ∀ t, x ∈ E c m : (cid:13)(cid:13) D m φ t π c m ( x ) (cid:13)(cid:13) ≤ c m , c e λ m , c | t | (cid:107) π c m ( x ) (cid:107) , (A3) ∀ t ≥ , x ∈ E s m : (cid:13)(cid:13) D m φ t π s m ( x ) (cid:13)(cid:13) ≤ c m , s e − λ m , s t (cid:107) π s m ( x ) (cid:107) , (A4) ∀ t ≤ , x ∈ E u m : (cid:13)(cid:13) D m φ t π u m ( x ) (cid:13)(cid:13) ≤ c m , u e λ m , u t (cid:107) π u m ( x ) (cid:107) . (A5)For the sake of simplicity, we restrict the constants λ m ,x ( x = c , s , u) to the absolute condition: λ x = sup m λ m ,x . Under this absolute condition, the condition 3 for φ t ∈ C k can be rewritten witha gap r ≥ − λ s < − rλ c ≤ ≤ rλ c < λ u . Then, M NHIM is C k for some k ≤ r , and the manifoldis called an eventually and absolutely k -NHIM. For this k -NHIM, there exists a vector field of˜ φ t and a manifold ˜ M NHIM that are C k –close, that is, the manifold remains a k -NHIM or persistsunder C k perturbation.The theorem was extended to non-compact NHIMs. In this case, Q is a bounded geometryand C k is accordingly changed into C k.xb,u . (See Ref. 68 for those definitions.) Appendix B: The NHIM Persistence Theorem for RDSs
For each d -dimensional stochastic path instance B t , one can uniquely obtain a d -dimensional, t -continuous function ω ( t ) —as defined in Appendix C)— that is associated with the saddle point ateach instance of time t , and for which we are free to initialize at ω (0) = 0. ω ( t ) is a generalizationof the so-called TS trajectory discussed in the main text. A probability distribution can thenbe defined for the bundle of instances ω from the probability distribution of B t , Thus, the event34pace Ω, which contains this bundle, is defined by the corresponding d -dimensional t -continuousfunctions C ( R , R d ): Ω = (cid:110) ω (cid:12)(cid:12)(cid:12) ω ∈ C ( R , R d ) , ω (0) = 0 (cid:111) . (B1)We define the Wiener shift θ s , and a cocycle φ t ω —that is a random invariant manifold M ( ω )—such that, φ t ω ( M ( ω )) = M ( θ t ω ). The persistence of the compact NHIM theorem for RDS holds for these invariant manifolds. Let us enumerate the differences in the expression of the persistencetheorems between Refs. 68 and 74.1. The following are all ω dependent: W (s)NHIM ( ω ), W (u)NHIM ( ω ), E x ( ω ), π x m ( ω ), c m ,x ( ω ), and λ m ,x ( ω ) ( x = s , u , c), as well as M NHIM ( ω ), and φ t ω .2. Eq. (A2) holds only for E x ( x = c , u). For E s , D m φ t ω ( E s m ( ω )) ⊂ E s D m φ t ω ( m ) ( θ t ω ) . (B2)3. (Persistence) M NHIM ( ω ) is a 1-normally hyperbolic random invariant manifold.However, M NHIM ( ω ) is still C k ( r ≥ k ) smooth when φ t ω is C k and − λ s ( ω ) < − rλ c ( ω ) < For the Langevin type SDEs, one can obtain its RDEs from the use of a stationary orbit .Suppose, a SDE with a Wiener process W t ∼ N (0 , t ),d X t = ( aX t + b ( X t ))d t + c d W t , (C1)is transformed by a stationary orbit η t , such that,d η t = aη t d t + c d W t , (C2)thus for t ≤ t ≤ t , d( X t − η t )/d t = ( ax + b ( X t )) , (C3) η t = e a ( t − t ) η t + c (cid:82) tt e − a ( s − t ) d W s e a ( t − t ) η t − c (cid:82) t t e − a ( s − t ) d W s , (C4)35here we use Itˆo’s lemma: d (cid:0) e − at η t (cid:1) = d (cid:0) e − at (cid:1) η t + e − at d η t + d (cid:0) e − at (cid:1) d η t = e − at ( − a d t + d η t ), andwe only wrote terms with order d t or less. When x = X − η such that x ( t ) = X t and ω ( t ) := η t with η = 0, one obtain a corresponding RDE,˙ x = ( ax + b ( x + θ t ω )) . (C5)For this equation, one have the corresponding RDS. The similar discussion can be made for thehigher dimensional, Langevin type (Eq. (C1)) systems. Here one can observe the TS trajectory η ‡ t , η ‡ t = − c (cid:82) t −∞ e − a ( s − t ) d W s ( a < c (cid:82) ∞ t e − a ( s − t ) d W s ( a > , (C6)is a special case of the stationary orbit Eq. (C2) by replacing t → −∞ , t → ∞ in Eq. (C4). Appendix D: Sampling the Reactivity Boundary In this subsection, we demonstrate how the algorithm, introduced in Sec. III and Fig. 6, worksin 1 D ketene with external force (Eq. (28)). Ideally, it is better to just locate the asymptotictrajectories, and this can be achieved by the perturbation theories. However, for irregular reac-tions, there could be a case for which the theories are not applicable, and need to use numericalinvestigation until an applicable theory is developed.The key step in the algorithm is the minimization of the uniform sampling. In Fig. D.2, wevisualize the locating process. Fig. D.2a, shows initial conditions by (2 + 1) × (2 + 1) gridsampling. This produce more than one sample for each area divided by the reactivity boundariesthus the sampling is sufficient as seeds for the algorithm. In the next step (Fig. D.2b), we apply thebisection method until when the demanded number of samples are obtained from the resolution.In the figure, the bisection method is only applied once for the visualization purpose. Then in thebottom figure, the results of the search on the resolution is given by the use of the algorithm weintroduced in Fig. 6. Finally, one can apply the bisection method for each pair until when thedemanded precision is achieved.Fig. D.3 shows the computational costs of the different sampling resolutions. That is, when wechose the sampling resolution (2 N res + 1) × (2 N res + 1) (Fig. D.2-bisection), and chose the precisionachieved by (2 + 1) × (2 + 1) grid, the x –axis of the figure is given by 2 N res + 1, and y –axis isgiven by actual computational time observed by std::clock , which is implemented in the C++standard template library. In the log-log plot, the cost is an almost linear-order increase over the36 IG. D.2. Step-by-step visualization of the algorithm to sample neighboring points of the forward-timereactivity boundaries in ketene 1 D with external force. The solution lines are shown in light gray. Straddlingpairs of the forward-time reactivity boundaries (gray links with dots). (a) The initial conditions absorbedat right (left) of each cell divided by the coordinate position of q F = − 3, min1, min2, min3, and q F = 3after positive time integration (blue (red) dots). (b) The initial straddling pairs before (after) the bisectionmethod applied (light gray (dark gray)). (c) The result of the path search algorithm (see Fig. 6) and itsmidpoint approximation of the neighboring points (blue). sampling resolution. The order is able to estimate by the parameter a ≈ y = a log x + b where ( a, b ) = (1 . , . ω of the periodic field, the phase space at t = 0 is identical when ωt ≡ π ). To take this advantage, one can propagate time t andget more samples. In Fig. D.4, we depict the results of time ( | ωt | = 0 , π, π, π ) integrations of37 number of samples along an axis10 − c o m pu t a t i o n c o s t s / s ec log y = a log x + b FIG. D.3. Computation cost in seconds as a function of the sampling resolution (see also Fig. D.2). Theinitial grid size is (2 + 1) × (2 + 1). The first bisection is applied up to (2 N res + 1) × (2 N res + 1) gridresolution where (2 N res + 1) is value of the x -axis. After the path search, the bisection method is appliedagain up to (2 + 1) × (2 + 1) grid resolution. The results are shown as cross symbols. The least squarefit is log y = 1 . × log x + 9 . the straddling pairs of manifolds obtained by the algorithm. In the figure, the propagated samplepoints (gray), midpoint estimation of each pair, and the line between the estimations (red or orangefor τ b , blue or cyan for τ f ) are shown. The strength of each color is proportional to the phase | ωt | .For example for the phase | ωt | = 2 πM , M = 0 is darkest, and the color get lighter for larger M . Weemploy N = 12 for the sampling resolution (2 N res + 1) × (2 N res + 1) of the algorithm. The gray dotsand lines are not visible because the estimation is precise, and those are overwritten by midpointestimations. The lines between the midpoints are given when the two consecutive midpoints haveless than a certain distance on the figure. The disconnection appears from the phase M ≥ e.g. for | ωt | ≡ π (mod 2 π ). In Fig. D.5,we show the manifolds at the phase | ωt | = π + 2 πM ( M = 0 , . . . , 6) (colored as the same manner asin Fig. D.4) superimposed on the phase | ωt | = 2 πM ( M = 0 , . . . , 7) (gray). Similarly, an arbitraryintermediate phase can be obtained. 38 IG. D.4. The reactivity boundaries at phase | ωt | = 0 , π, π, π (a-d) obtained by a method explainedin Fig. D.2 with the resolution given by N res = 16. The reactivity boundaries are corresponding to stable(blue, cyan) and unstable (red, orange) manifolds. The colors of the manifolds get lighter as | t | increases. ∗ [email protected] H. Eyring, J. Chem. Phys. , 107 (1935), doi:10.1063/1.1749604. M. G. Evans and M. Polanyi, Trans. Faraday Soc. , 875 (1935), doi:10.1039/tf9353100875. E. P. Wigner, Trans. Faraday Soc. , 29 (1938), doi:10.1039/TF9383400029. S. Glasstone, K. J. Laidler, and H. Eyring, The theory of rate processes : the kinetics of chemicalreactions, viscosity, diffusion and electrochemical phenomena (McGraw-Hill Book Company, inc., NewYork and London, 1941). IG. D.5. The stable (unstable) manifolds at phase | ωt | ≡ π (mod 2 π ) (blue, cyan (red, orange)). Acomparison can be made with another phase | ωt | ≡ π ) (gray) in the background. K. Fukui, J. Phys. Chem. , 4161 (1970), doi:10.1021/j100717a029. S. Kato and K. Fukui, J. Am. Chem. Soc. , 6395 (1976), doi:https://doi.org/10.1021/ja00436a061. K. J. Laidler and M. C. King, J. Phys. Chem. , 2657 (1983), doi:10.1021/j100238a002. D. G. Truhlar, W. L. Hase, and J. T. Hynes, J. Phys. Chem. , 2664 (1983), doi:10.1021/j100238a003. P. H¨anggi, P. Talkner, and M. Borkovec, Rev. Mod. Phys. , 251 (1990), and references therein,doi:10.1103/RevModPhys.62.251. D. G. Truhlar, B. C. Garrett, and S. J. Klippenstein, J. Phys. Chem. , 12771 (1996),doi:10.1021/jp953748q. J. C. Keck, Adv. Chem. Phys. , 85 (1967), doi:10.1002/9780470140154.ch5. F. T. Wall, L. A. Hiller, and J. Mazur, J. Chem. Phys. , 255 (1958), doi:10.1063/1.1744471. J. S. Wright, J. Chem. Phys. , 720 (1978), doi:10.1063/1.436639. E. Pollak and P. Pechukas, J. Chem. Phys. , 1218 (1978), doi:10.1063/1.436658. Y. Nagahata, H. Teramoto, C.-B. Li, S. Kawai, and T. Komatsuzaki, Phys. Rev. E , 042923 (2013),doi:10.1103/PhysRevE.88.042923. S. Patra and S. Keshavamurthy, Phys. Chem. Chem. Phys. , 4970 (2018), doi:10.1039/C7CP05912D. C. Dellago, P. Bolhuis, F. S. Csajka, and D. Chandler, J. Chem. Phys. , 1964 (1998),http://gold.cchem.berkeley.edu/Pubs/DC150.pdf. C. Dellago, P. G. Bolhuis, and P. L. Geissler, Adv. Chem. Phys. , 1 (2002),doi:10.1002/0471231509.ch1. T. Bartsch, R. Hernandez, and T. Uzer, Phys. Rev. Lett. , 058301 (2005),doi:10.1103/PhysRevLett.95.058301. R. Hernandez, T. Bartsch, and T. Uzer, Chem. Phys. , 270 (2010),doi:10.1016/j.chemphys.2010.01.016. M. J. Davis, Chem. Phys. Lett. , 491 (1984), doi:10.1016/0009-2614(84)87077-3. R. Mackay, J. Meiss, and I. Percival, Physica D , 55 (1984), doi:10.1016/0167-2789(84)90270-7. A. M. Ozorio de Almeida, N. de Leon, M. A. Mehta, and C. C. Marston, Physica D , 265 (1990), oi:10.1016/0167-2789(90)90040-V. S. Wiggins, Physica D , 471 (1990), doi:10.1016/0167-2789(90)90159-M. S. Wiggins, L. Wiesenfeld, C. Jaffe, and T. Uzer, Phys. Rev. Lett. , 5478 (2001),doi:10.1103/PhysRevLett.86.5478. N. Fenichel, Indiana Univ. Math. J. , 193 (1972), doi:10.1512/iumj.1972.21.21017. R. Hernandez and W. H. Miller, Chem. Phys. Lett. , 129 (1993), doi:10.1016/0009-2614(93)90071-8. R. Hernandez, J. Chem. Phys. , 9534 (1994), doi:10.1063/1.467985. T. Komatsuzaki and M. Nagaoka, J. Chem. Phys. , 10838 (1996), doi:10.1063/1.472892. T. Uzer, C. Jaff´e, J. Palaci´an, P. Yanguas, and S. Wiggins, Nonlinearity , 957 (2002), doi:10.1088/0951-7715/15/4/301. H. Waalkens, R. Schubert, and S. Wiggins, Nonlinearity , R1 (2008), doi:10.1088/0951-7715/21/1/R01. S. Kawai and T. Komatsuzaki, J. Chem. Phys. , 084304 (2011), doi:10.1063/1.3554906. ¨U. C¸ ift¸ci and H. Waalkens, Nonlinearity , 791 (2012), doi:10.1088/0951-7715/25/3/791. S. Kawai and T. Komatsuzaki, J. Chem. Phys. , 224505 (2009). S. Kawai and T. Komatsuzaki, J. Chem. Phys. , 224506 (2009). S. Kawai and T. Komatsuzaki, Phys. Chem. Chem. Phys. , 15382 (2010), doi:10.1039/c0cp00543f. S. Kawai, A. D. Bandrauk, C. Jaff´e, T. Bartsch, J. Palaci´an, and T. Uzer, J. Chem. Phys. , 164306(2007), doi:10.1063/1.2720841. S. Kawai and T. Komatsuzaki, J. Chem. Phys. , 024317 (2011), doi:10.1063/1.3528937. C.-B. Li, A. Shoujiguchi, M. Toda, and T. Komatsuzaki, Phys. Rev. Lett. , 028302 (2006),doi:10.1103/PhysRevLett.97.028302. S. Kawai and T. Komatsuzaki, Phys. Rev. Lett. , 048304 (2010),doi:10.1103/PhysRevLett.105.048304. T. Bartsch, T. Uzer, J. M. Moix, and R. Hernandez, J. Chem. Phys. , 244310 (2006),doi:10.1063/1.2206587. D. Townsend, S. A. Lahankar, S. K. Lee, S. D. Chambreau, A. G. Suits, X. Zhang, J. L. Rheinecker,L. B. Harding, and J. M. Bowman, Science , 1158 (2004), doi:10.1126/science.1104386. J. M. Bowman and B. C. Shepler, Annu. Rev. Phys. Chem. , 531 (2011), doi:10.1146/annurev-physchem-032210-103518. F. A. Maugui`ere, P. Collins, Z. C. Kramer, B. K. Carpenter, G. S. Ezra, S. C. Farantos, and S. Wiggins,Annu. Rev. Phys. Chem. , 499 (2017), doi:10.1146/annurev-physchem-052516-050613. I. S. Ulusoy, J. F. Stanton, and R. Hernandez, J. Phys. Chem. A , 7553 (2013),doi:10.1021/jp402322h. H. Teramoto, M. Toda, and T. Komatsuzaki, Phys. Rev. Lett. , 054101 (2011),doi:10.1103/PhysRevLett.106.054101. Y. Nagahata, H. Teramoto, C.-B. Li, S. Kawai, and T. Komatsuzaki, Phys. Rev. E , 062817 (2013), oi:10.1103/PhysRevE.87.062817. G. Haller, Annu. Rev. Fluid Mech. , 137 (2015), doi:10.1146/annurev-fluid-010313-141322. A. Hadjighasem, M. Farazmand, D. Blazevski, G. Froyland, and G. Haller, Chaos , 053104 (2017),doi:10.1063/1.4982720. C. Mendoza and A. M. Mancho, Phys. Rev. Lett. , 038501 (2010),doi:10.1103/PhysRevLett.105.038501. G. T. Craven and R. Hernandez, Phys. Rev. Lett. , 148301 (2015),doi:10.1103/PhysRevLett.115.148301. A. Junginger and R. Hernandez, J. Phys. Chem. B , 1720 (2016), doi:10.1021/acs.jpcb.5b09003. G. T. Craven and R. Hernandez, Phys. Chem. Chem. Phys. , 4008 (2016), doi:10.1039/c5cp06624g. F. Revuelta, R. M. Benito, and F. Borondo, Phys. Rev. E , 032221 (2019),doi:10.1103/PhysRevE.99.032221. J. D. Gezelter and W. H. Miller, J. Chem. Phys. , 7868 (1995), doi:10.1063/1.470204. K. Ide, D. Small, and S. Wiggins, Nonlinear Process. Geophys. , 237 (2002), doi:10.5194/npg-9-237-2002. J. A. Jim´enez Madrid and A. M. Mancho, Chaos , 013111 (2009), doi:10.1063/1.3056050. A. M. Mancho, S. Wiggins, J. Curbelo, and C. Mendoza, Commun. Nonlinear Sci. Numer. Simul. ,3530 (2013), doi:10.1016/j.cnsns.2013.05.002. G. Haller, “Interactive comment on ‘Detecting and tracking eddies in oceanic flow fields: A vorticitybased Euler-Lagrangian method’ by R. Vortmeyer-Kley et al.” (2016), doi:10.5194/npg-2016-16-SC2. F. T. Wall, L. A. Hiller, and J. Mazur, J. Chem. Phys. , 1284 (1961), doi:10.1063/1.1732040. F. T. Wall and R. N. Porter, J. Chem. Phys. , 3112 (1963), doi:10.1063/1.1734151. J. S. Wright, G. Tan, K. J. Laidler, and J. E. Hulse, Chem. Phys. Lett. , 200 (1975), doi:10.1016/0009-2614(75)80100-X. J. S. Wright, K. G. Tan, and K. J. Laidler, J. Chem. Phys. , 970 (1976), doi:10.1063/1.432291. J. S. Wright and K. G. Tan, J. Chem. Phys. , 104 (1977), doi:10.1063/1.433656. K. J. Laidler, K. Tan, and J. S. Wright, Chem. Phys. Lett. , 56 (1977), doi:10.1016/0009-2614(77)85162-2. K. G. Tan, K. J. Laidler, and J. S. Wright, J. Chem. Phys. , 5883 (1977), doi:10.1063/1.434795. S. Wiggins, Normally Hyperbolic Invariant Manifolds in Dynamical Systems (Springer New York, NewYork, NY, 1994) doi:10.1007/978-1-4612-4312-0. J. Eldering, Normally Hyperbolic Invariant Manifolds (Atlantis Press, Paris, 2013) doi:10.2991/978-94-6239-003-4. V. I. Arnol’d, Mathematical Methods of Classical Mechanics , 2nd ed., edited by J. H. Ewing, F. W.Gehring, and P. R. Halmos, Graduate Texts in Mathematics, Vol. 60 (Springer-Verlag, UNKNOWN,1989). M. W. Hirsch, C. C. Pugh, and M. Shub, Invariant Manifolds (Springer Berlin Heidelberg, Berlin, eidelberg, 1977) doi:10.1007/BFb0092042. B. Øksendal, Stochastic Differential Equations: An Introduction with Applications (Springer Berlin Hei-delberg, Berlin, Heidelberg, 2003) doi:10.1007/978-3-642-14394-6. J. Duan, An Introduction to Stochastic Dynamics , 1st ed., Cambridge Texts in Applied Mathematics(Cambridge University Press, 2015) p. 312. P. E. Kloeden and E. Platen, Numerical Solution of Stochastic Differential Equations (Springer BerlinHeidelberg, Berlin, Heidelberg, 1992) doi:10.1007/978-3-662-12616-5. J. Li, K. Lu, and P. Bates, Transactions of the American Mathematical Society , 5933 (2013),doi:10.1090/S0002-9947-2013-05825-4. L. Arnold, Random Dynamical Systems , 2nd ed. (Springer, Berlin, Heidelberg, 2003) doi:10.1007/978-3-662-12878-7. I. S. Ulusoy and R. Hernandez, Theor. Chem. Acc. , 1528 (2014), doi:10.1007/s00214-014-1528-z. A. P. Scott, R. H. Nobes, H. F. Schaefer III, and L. Radom, J. Am. Chem. Soc. , 10159 (1994),doi:10.1021/ja00101a039. E. R. Lovejoy, S. K. Kim, R. A. Alvarez, and C. B. Moore, J. Chem. Phys. , 4081 (1991),doi:10.1063/1.460764. E. R. Lovejoy, S. K. Kim, and C. B. Moore, Science , 1541 (1992), doi:10.1126/science.256.5063.1541. E. R. Lovejoy and C. B. Moore, J. Chem. Phys. , 7846 (1993), doi:10.1063/1.464592. “C++ boost libraries,” , accessed: 2019-01-30., accessed: 2019-01-30.