Phase space geometry of isolated to condensed chemical reactions

Questions about how and why chemical species change remain critical to all molecular sciences. The mechanism involves the reorganization of electrons and the motion of nuclei along a pathway typically discussed in their configuration space.¹ The theory of reaction rates dates back to van 't Hoff² and Arrhenius³ and was elaborated by Marcelin.^4,5 It was later reformulated as transition state theory^6–9 and Rice–Ramsperger–Kassel–Marcus (RRKM) theory.^10–14 Advances in classical,^15–27 semiclassical,^28–33 and quantum mechanical³⁴ reaction rate theories have aimed to recast this approximate (but generally accessible) theory to obtain accurate, and sometimes exact, rates. Formally, the requisite quantum mechanical formulas to obtain correct and direct rates are known as long as the equations of motion can be integrated correctly. Unfortunately, such exact quantum mechanical numerics are not yet within reach for the full all-atom treatment of chemical reactions in solvents.

A key concept introduced in all these rate theories is the transition state (TS) that acts as a doorway between reactants and products. Thereby, the process of reaction is decomposed into two subprocesses: when the system wanders about the reactant (or product) before finding the TS and when the system crosses the TS after having first reached it. The former process has often been assumed to be locally “equilibrated” before the system finds the TS. This may be understood intuitively for isolated systems. If the total number of degrees of freedom (DoF) n is much greater than unity (n ≫ 1) and the total energy E is held fixed, then the likelihood for the system to find the TS may be rare simply because of the search problem. In addition, the barrier height E_saddle (≫E/n) is generally large compared to the average energy available to the reactive motion, and this further reduces the likelihood that the system will access the TS. In turn, within the TST framework, the latter process has been assumed such that once the system finds and crosses the TS, it will not return to the TS before being “captured” (viz, “equilibrated”) in the product. However, the accuracy of this approximation is questionable because for most naive TSs, the correct rates can only be captured in the rate formula by way of introducing an ad hoc correction factor, known as the transmission coefficient κ(≲1). Since the early work of Eyring^6,8 and Wigner,^7,9 this factor has been known to account for nonlinear and dynamical effects (in addition to quantum, tunneling effect) not being easily accessible to statistical thermodynamic considerations. For example, κ takes into account the existence of recrossing trajectories that generally lead the transition state theory (TST) to overestimate the reaction rate constant.

Nevertheless, the underlying beauty of TST is that it replaces the formally complex and numerically expensive calculation of trajectories with the introduction of a geometric object—the TS—in between reactants and products. In some past treatments of chemical reaction theory, the TS has also been regarded as a quasi-stable “activation complex” lying somewhere between the reactant and the product in the configuration space. In this Perspective, the TS is defined more narrowly as a dividing surface (DS) from which one can obtain a flux from the reactant states to the product states. This definition is rooted in the early work of Wigner⁹ and Keck¹⁵ following Marcelin.^4,5 The most striking consequence of the phase space perspective overlaid onto the theory of chemical reactions is not just that it leads to accurate rates of reactions (in terms of the exact TS), but also that it helps to characterize the mechanism of why and how a reaction takes place.

A spacecraft moving through one planetto another can also be regarded as a “reaction” although its spatial scale is different (at ∼10²⁰ times larger than that of chemical reactions). The stable states characteristic of a chemical reaction correspond to trapped regions due to the net gravitational interactions. For example, the escape rate of an asteroid from a planet to another planet is a dynamical transition that can be characterized quantitatively using the concept of a TS.^35,36 Moreover, in 2000,³⁶ pathways for spacecraft to move between planets were found to exist under the influence of their gravitational forces, just as a “hang glider” moves through the air. This greatly expanded the possible route designs in space beyond “brute force” driven only by fuel and has the potential to drastically reduce the amount of fuel and weight of spacecraft.³⁷ A critical structure to uncovering this pathway is the stable and unstable invariant manifolds emanating from the normally hyperbolic invariant manifold (NHIM),^38–42 which form the boundary of reactivity to precisely divide trajectories into reactive and nonreactive ones (viz, reactivity boundary), from dynamical systems theory. In the space mission problem, a spacecraft moves in the space where planets are approximated as “being frozen” in the rotational coordinate system so that the effective DoF is maximally three if the spacecraft is represented as a point. In contrast, chemical reactions are composed of atoms (whose masses are comparable) and are essentially multidimensional, which provides several challenging subjects not only as the problems in chemistry but also as those in dynamical systems theory: e.g., the question of the bifurcation and/or breakdown of NHIM^{19–23,43–50} in reference to the origin of a stochastic process, i.e., when and how the phase space geometry may be “broken” so that chemical reactions essentially stochastic and unpredictable.

Not surprisingly, chemists have explored the phase space properties of the TS in the context of chemical reactions. In the limit of two DoF, such reactions can be represented through the potential energy surface as shown in Fig. 1. Approaches for resolving the dynamics with respect to geometrical structures associated with the reactant, product, and in-between regions include the reactivity band/map theory^51–56 (initiated by Wall et al.⁵⁷), the periodic orbit dividing surface (PODS) corresponding to the projection of NHIM onto the coordinate space,^19–23 reactive island theory,^25–27 and turnstile theory.^58,59 Unfortunately, despite a few notable exceptions—e.g., Refs. 60 and 61, which addressed the possibility of the nonexistence of turnstiles in three and higher DoF—these efforts were often limited to two-DoF systems due to its high computation cost. Thus, these phase space-based theories grew out of favor as brute-force higher-dimensional molecular dynamics calculations became accessible. In the last two decades, however, the use of transition state theory leveraging phase space perspectives have received renewed interest in the chemical literature. The existence of a reaction pathway through which all reactive trajectories necessarily follow, just like in the space mission route design,³⁶ has been revealed in a wide class of chemical reactions of many-DoF systems in the gas^{43,44,62–94} and condensed phases.^95–111 They also gave rise to advances in the corresponding semi-classical theories,^31,112–114 which provides a solid interpretation to the question of why an optimally controlled laser field designed to make a reaction is so complicated¹¹⁴ (see reviews in Refs. 115–123).

In this article, we review our current understanding of the generalization of the PODS in relation to the NHIM—also illustrated in Fig. 1—beyond our earlier reviews^{115,119–123} and its potential applications to gas and condensed phase chemical reactions. These approaches are complementary to those that emphasize the analysis of a reduced set of trajectories, such as transition paths,^124–129 in that such trajectories must somehow cross the DS. The purpose of this Perspective is thus to accelerate the bridging between the mathematical and chemical disciplines in using phase space geometric structures to more accurately determine rates and mechanisms in high-dimensional chemical reactions.

The phase space geometry of generic two-DoF Hamiltonian systems in the vicinity of an index-one saddle is illustrated in Fig. 2. Here, $(p1′,q1′)$ and $(p2′,q2′)$ are reactive and nonreactive normal form coordinates and conjugate momenta, respectively.¹³⁰ The nonlinear couplings between the normal modes can be locally eliminated by iterating perturbative corrections. When the total energy E of the system is just slightly above the potential energy at the saddle point E_saddle, this procedure does not destroy the normal mode structure though it is clearly transformed. Thus, $(p1′,q1′)$ and $(p2′,q2′)$ can be regarded as being independent of each other in the vicinity of the saddle, giving rise to an unstable periodic orbit (UPO). The UPO in phase space corresponds to a closed loop [shown as a thick purple curve denoted by $M$ in Fig. 2(c)]. A bundle of trajectories asymptotically approaching this closed loop but never reaching it while remaining in its proximity and those leaving the loop asymptotically at infinite time [shown as a filled purple circle and lines in Figs. 2(a) and 2(b)] together define a structure of interest. For an index-one saddle of an energy surface (in arbitrary dimensions), it has the topology of a tube in phase space. For more complex reaction geometries, the structure still exists locally as a boundary between reactive and nonreactive regions (known as the reactivity boundary), which we refer to as tube-like. When the system resides either on the closed loop or on the reactivity boundary (or tube in special cases) itself, it does so perpetually. In this sense, these phase space objects, the UPO and the reactivity boundary, serve as (time) invariant manifolds.

The reactive (red) trajectories in Fig. 2(a) have larger energies in the reactive mode than those of the reactivity boundary (purple lines), which is a tube in the case shown there. Similarly, the energy component associated with the projection of the reactive trajectories to $(p2′,q2′)$ should be smaller than those of the tube [purple circle in Fig. 2(b)] due to the conservation of the total energy. Hence, all trajectories that will react pass through the inside of the corresponding tube either from the reactant to the product or vice versa. No trajectory can cross the reaction tube because the tube forms an invariant manifold of the equation of the motion. In other words, the tube serves as “an impenetrable barrier” in phase space and does act as a reactivity boundary (as claimed above) between regions of reactive and nonreactive trajectories. It is also sometimes called a reactive island^25–27 [Fig. 2(c)].

The PODS was an important advance in the identification of non-recrossing DSs in the configuration space, which is defined as the projection of the UPO in the phase space onto the configuration space.^19–23 It was examined for some systems such as collinear triatomic exchange reactions of HX + H (with X = H, Cl, F) (a two-DoF system) that the projection of the UPO onto the configuration space forms non-recrossing DSs with the required dimensionality to divide that configuration space into two disjoint regions, i.e., reactant and product. The required dimensionality of DS is that any surface to divide the given space (where reaction events are represented) into those before and after reaction events should have one dimension less than that of the space. For example, in two (three)-dimensional space, the required dimensionality of the DS is the one (two)-dimensional manifold (one- and two-dimensional manifolds correspond to simply a line and a two-dimensional surface, respectively); otherwise, no one can see trajectories “crossing” the DS generically [e.g., a point (zero-dimensional manifold) does not divide the space of more than one-dimensional space]. It should be noted for n-DoF systems that the projection of a set of UPOs or more in general NHIM that are geometrical objects in 2n-dimensional phase space onto the n-dimensional configuration space has no guarantee to satisfy the required dimensionality with no-return property (see also some relevant studies^131,132).

In what remains of this section, we discuss the perturbative and nonperturbative approaches that have emerged to address these limitations of a no-return DS in the configuration space by lifting it into phase space.

Since the 1990s, a non-recrossing dividing surface has been extracted in phase space for many-DoF Hamiltonian systems using perturbative approaches.^{31,43,44,63–70,72–80,83} In order to preserve the Hamiltonian structure in phase space, the selected perturbation theory must be symplectic. The crux to doing so is the use of Lie canonical perturbation theory (Lie-CPT)^{130,133–139} for Hamiltonian systems H. It expands as a Taylor series in normal coordinates q around index-one saddle points in which the potential is concave down along only one direction, say, q₁,

where

Here, q_i and p_i are the ith normal coordinate and its conjugate momentum, respectively; ω_i and C_jkl, C_jklm, … are, respectively, the frequency of the ith mode and the coupling coefficient among q_j, q_k, and q_l and that among q_j, q_k, q_l, and q_m and so forth. The frequency associated with q₁, ω₁ ≡ − |ω₁|i (because of the negative curvature $ω12<0$ along q₁) and those of the other modes are pure-imaginary and real, respectively.

When the Hamiltonian H is the nonperturbed Hamiltonian, i.e., H = H₀,

or

Equation (4) can be rewritten symmetrically with respect to action-angle canonical variables (J, θ) as

resulting in the time-invariance of half of the canonical variables—viz, the action variables. As a consequence, the system is integrable, leading to analytical solutions of the equation of motion. It should be noted that in this formulation, the normal mode “action” and the conjugate angle variable along the hyperbolic reactive DoF 1 are not given by the standard definition $Ji≔12π∮pidqi(=12ωi(pi2+ωi2qi2))$ and Θ_i =  tan⁻¹(ω_iq_i/p_i). The action and angle variables associated with the reactive DoF is formulated by $J1≔i2|ω1|(pi2−|ω1|2qi2)$ and $Θ1=−i⁡tanh−1|ω1|q1p1$⁠. It is easily verified (e.g., Ref. 70) that those action and angle variables including those of nonreactive DoF are canonical although there exists no period along hyperbolic mode. This is an adoption of Miller’s action variable²⁸ for the estimation of the tunneling reaction: $J1≔∫ab2(V1(q1)−E1)dq1$⁠, where V₁ and E₁ are the potential energy and mode energy associated with q₁, respectively, and a, b are points satisfying V₁(q₁) = E₁.

CPT canonically transforms a coordinate system (p, q) into a new coordinate system (p′, q′) so as to make the Hamiltonian simpler (e.g., decoupled among new modes) as much as possible. The transformation is subsequently operated as order by order so that the new Hamiltonian is made successively more integrable to increasing order. Specifically, when CPT can transform a given Hamiltonian H(p, q) into a new Hamiltonian $H′=H′J′+O(ϵk)$⁠, where J′ is the new action variables associated with (p′, q′) so that $dJi′/dt=Oϵk$ (i.e., $Ji′$ is approximately time invariant), the equations of motion in the old phase space coordinates are transformed^64,70 into

Thus,

where the transformed new frequency $ωi′(J′)$ is given by time derivative of the new angle variable $Θi′$ conjugate to the new action variable $Ji′$⁠: $ωi′(J′)=dΘi′/dt=∂H′(J′)∂Ji′=const+Oϵk$⁠. The transformed variables $(pi′,qi′)$ are approximately independent of each other as long as $ωi′$ is regarded as constant in time (cf. $d2qi/dt2=−ωi2qi$ for the normal mode Hamiltonian system where each mode is decoupled) as illustrated in Fig. 3. One can derive analytical expressions of (p′, q′) with respect to (p, q) and vice versa (see the full derivation in Ref. 70).

CPT is regarded as a normal form theory^{130,140–142} for Hamiltonian systems and a classical analog of Van Vleck perturbation theory. It was applied to obtain good-action variables in semi-classical TST by Miller and his co-workers.^{28–33,140,141,143,144} Hernandez and Miller pointed out³² that resonances that can spoil such perturbation theories are never met among nonreactive modes and the reactive mode because their associated frequencies are real and imaginary, respectively. In Refs. 31 and 145, they further observed that this construction of semiclassical TST could be used to obtain improved DSs in the classical cases of TST discussed here. Komatsuzaki et al.^63,64 showed by applying Lie-CPT with algebraic quantization to a four DoF Hamiltonian system of intramolecular proton transfer reaction of malonaldehyde that TST generalized in terms of no-return TS $q1′(p,q)=0$ can analytically provide the ad hoc transmission coefficient κ that used to be determined empirically. They^65–69,115 derived the analytical expression of a reaction coordinate $q1′(p,q)$ decoupled from the rest under the existence of chaos in a nonreactive space¹¹⁵ following Ref. 64 and showed that a no-return TS can persist even though the nonreactive coordinates cannot be decoupled in a six argon atom cluster by the second order Lie-CPT. This implies that the impossibility of a resonance among nonreactive modes and the reactive mode³² classically corresponds to the case that new nonreactive modes cannot be decoupled with each other due to the resonance among those modes, but they can be still effectively decoupled from the new reactive mode $q1′$⁠. They also addressed the regularity of saddle crossing dynamics dependent on the total energy above the saddle point energy, such as quasi-regular, semi-chaotic, and fully-developed chaos regions. Wiggins and his co-workers^71,72,112 provided the underlying phase space geometrical structure of chemical reactions in terms of the NHIM, corresponding to the generalization of reactive island theory^25–27 to many-DoF systems, and addressed the importance of dimensionality required for TS to divide the ambient space into two disjoint subspaces of reactant and product and the nonuniqueness of the no-return TS.

We obtain an analytical formulation of the DS as

not in the configuration space but in the phase space, which is free from recrossings because there exists no coupling among those modes at least up to $Oϵk$ in the vicinity of saddle. The analytical expression of the reaction coordinate $q1′(p,q)$ implies that one can predict the fate of reactions (i.e., the trajectory in question will react or not) without calculating trajectories once we know (p, q) at some transit time.

When canonical normal form is estimated, the geometry of the structures described above can be captured through the following building blocks: the NHIM, $M$⁠,

and its stable and unstable invariant manifolds, $Ws$ and $Wu$⁠, expressed as

and

respectively. In two-DoF systems, the NHIM simply corresponds to the UPO and the stable and unstable invariant manifolds $Ws$ and $Wu$ to the boundary of the reaction tube shown in Fig. 2. The underlying rigorous mathematics of the phase space geometry of the NHIM and its associated stable/unstable invariant manifolds are reviewed in Refs. 71, 72, and 112 and for chemical physics community in Refs. 115–123. The NHIM can be regarded as a “generalized saddle” in the phase space of many-DoF systems because the system perpetually stays in the NHIM while the system has non-zero momenta.

The NHIM and its stable and unstable invariant manifolds have an important property: they are persistent under sufficiently small and smooth perturbations, such as from moderate changes in the system’s total energy, effects from an external field, or deviations in the potential arising from the level of approximation in electronic structure calculations. In other words, these invariant manifolds change only smoothly—i.e., do not disappear suddenly—under such small perturbations.^38–42 This means that if one manages to obtain a transformation to some order, say, k − 1, then the NHIM constructed by that transformation gives the corresponding manifold because it persists under the map of a k-differentiable function. This property (called “structurally stability”) has not received much attention in the chemical physics literature but it is useful here. Mathematically, it guarantees that even while we are invariably dealing with an approximate Hamiltonian in describing a chemical system (e.g., due to the level of approximation of electronic structure theory or limits on the accuracy of the force field one may use), we can still obtain qualitatively accurate information about the true dynamics (which result from the propagation of the exact Hamiltonian). The other important property concerns the dimensionality of these manifolds. For instance, the phase space TS defined by Eq. (9) is of codimension one, which is one dimension less than that of the full space. This property is required to divide a given space into two disjoint subspaces, such as the regions of the reactant or product; e.g., in order to divide a three (two)-dimensional space into two regimes, one needs a two-dimensional plane (one-dimensional line). Following Wigner’s definition of the TS, it is a surface at the “point of no-return” through which all reactive trajectories cross only once. In order for the term “cross” to make sense, the TS should be of co-dimension one; otherwise, one cannot count the event of crossing in the space as that when the reaction occurs.

The NHIM and the reaction tube itself are invariant and uniquely defined at a given total energy of the system E, but the TS, defined as a co-dimension one surface free from recrossing, is not unique. Such multiplicity arises because any “disk” or “valve” within the reaction tube can serve as a TS (e.g., yellow line in Fig. 2). When the total energy of a system E is not far above from the saddle point energy E_saddle and the nonlinear couplings among the modes can be treated as small perturbations, one can generalize the concept of the no-return DS to many-DoF Hamiltonian systems even in the chaotic regime. However, when the total energy of a reacting system increases beyond this, this construction is no longer so simple.^{19–23,43–50,83}

First, as seen in the pioneering work developing the PODS,^19–23 the NHIM (and hence, the associated TSs) can bifurcate as a function of the total energy of the system in many-DoF systems. The challenging complexity beyond the case of two-DoF systems is that the control parameters of bifurcation are no longer just the total energy but also the action variables. If and only if the action variables are approximately preserved near the regions where the NHIM bifurcates, then one can draw the bifurcation diagram of the NHIM as a function of the bath modes normal form action(s) and predict the bifurcation of the NHIM.⁴⁴ Second, as shown for a hydrogen atom in crossed electric and magnetic fields,^46,49,50 as the energy increases further, a no-return TS may no longer exist, but the NHIM and the stable (unstable) invariant manifolds that serve as “an impenetrable barrier” in the phase space persist, just like the emergence of “turbulence” within a reaction tube does not necessarily destroy the tube.⁸³ 

Other crucial directions along the perspective afforded by the phase space geometry of chemical reactions are as follows: How does the phase space geometry in classical mechanics affect the corresponding quantum mechanics? The answers to this question may still be far from being fully satisfied, but the corresponding semi-classical theories have been developed^112–114 following the pioneering works.^146–148 For example, it was shown that optimally controlled,^149,150 complicated laser field that brings a vibrational state localized on one potential minimum to another vibrational state localized on another potential minimum was simply interpreted as the “catch-and-release mechanism” by the underlying reaction boundary modulated by the time-dependent laser field.^113,114 How can one incorporate the spatial rotation coupled with the internal vibrational motions, and what type of new feature may arise in defining no-return TS and the reaction tube? The corresponding Hamiltonian is no longer as simple as Eqs. (1)–(3), but the mathematical formulation was made under the existence of rovibrational couplings.^151–153 Crossing dynamics through higher rank saddles, which is expected to be more important for phase transition-like motion in atomic clusters¹⁵⁴ and/or at higher energies compared to the saddle point energy, was also studied in terms of phase space geometry^85–89 but not explicitly using the NHIM.^88,89 In order to include this structure, one also needs to specify the role of transition states in determining the dynamics in a roaming reaction,^155–157 a valley-ridge bifurcation of reaction path,^158,159 the dynamical switching of a reaction coordinate,^46,49,50 and in an emerged transition state via external field.^160,161 The perturbation theories can provide the (approximate) analytic formula of the phase space reaction coordinate locally decoupled to the rest, yielding (approximate) no-return TS and stable and unstable manifolds (reaction boundary). However, those require a zeroth order approximation of a perturbation theory or are approximations acquired up to some order $Oϵk$⁠. There are several cases when the CPT approach faces challenges in convergence or otherwise, for example, when the total energy increases so much that nonlinear couplings among original coordinates may no longer be regarded as a perturbation. Another case arises when an increase in the dimensionality of the system leads the number of monomial terms to grow so rapidly (in a combinatorial manner) that the computation is infeasible.

Just recently, nonperturbative methods exploring the larger dimensional phase space were also developed in the context of reactions.^{46,49,50,110,160,162–164} For example, a hydrogen atom in crossed electric and magnetic fields was found to have a significant instability in its dynamics within the NHIM at the total energy of the system so high that a perturbative treatment is impossible.^46,49,50 As the total energy increases, the hyperbolicity along the reaction coordinate in the phase space starts to compete with the hyperbolicity developed within the NHIM. The normal hyperbolicity then begins to break down, and the NHIM can no longer persist though some remnants of its structure remain visible. As the total energy increases further, the hyperbolicity developed within the “NHIM” starts to dominate the hyperbolicity along the coordinate that used to be the reaction coordinate at the lower energy. A new normally hyperbolic structure emerges in which the old reaction coordinate has switched to a non-reaction coordinate. In Subsection II C, we now focus on some of the nonperturbative approaches that have been pursued as possible alternatives.

As summarized thus far, the powerful tools of dynamical systems theory have been used to resolve the geometry of chemical reactions since the late 1980s.^25,58,60 In many cases, the geometry—represented through normal forms—of the reaction could be resolved using perturbative methods. However, perturbative approaches become inapplicable when perturbation terms get larger and/or DoF of the system in question increases much more than, say, 100. This has led to the development of several alternative nonperturbative approaches for constructing the NHIM and related structures, such as the stable/unstable invariant manifolds.

One set of nonperturbative approaches for identifying the NHIM for chemical reactions has centered on the use of the so-called Lagrangian Descriptor (LD). A key associated concept is the Lagrangian coherent structure (LCS),^165,166 which is an invariant manifold in the Lagrangian frame—that is, in the moving frame associated with the phase space transport—and thereby provides a coordinate system tied to trajectories evolving in time. The LD^167–169 has been seen to be useful for uncovering the LCS.¹⁶⁶ The LDs first introduced in Ref. 167 heuristically support the conjecture that stable/unstable manifolds can be described by the “abrupt change”¹⁶⁸ of $M$-function (e.g., arc length of a trajectory over τ time evolution). Mendoza and Mancho¹⁶⁷ found that in a system with one DoF hyperbolic normal mode, a typical LCS descriptor—that is, finite time Lyapunov exponent ridge—was not able to uncover the stable and unstable manifold found using the LD.¹⁶⁸ Despite the controversy around the conjecture due to alleged counterexamples^170,171 and their possible resolution¹⁶⁹ and the fact that the LD is not objective¹⁷² (that is, the LD does not have invariance for a coordinate transformation), the LD has been seen to be numerically effective in an increasing number of dynamical problems.

In the context of reaction dynamics, the LD was introduced¹¹⁰ as the extremal LD $M$⁠, describing the arc length of trajectories in coordinate space,

where the forward (backward) LD $Lf$ (⁠$Lb$⁠) is defined as follows:

where t₀ is the initial time, q (⁠$q̇$⁠) is the coordinate (velocity) vector of the system, q₀ (⁠$q̇0$⁠) is that of the initial condition at time t₀, and ‖x‖_p is the L^p norm of vector x defined by $‖x‖p≔∑ixip1/p$ (L² is the Euclidean distance) of the velocity. As discussed in Ref. 168, $Lf$ and $Lb$ are the accumulation of a positive scalar along a given trajectory. χ^(f) and χ^(b) are renormalization factors accounting for the differences in timescales. For example, the characteristic exponents (expansion and contraction ratio) of the Langevin equation (LE) are different and can be corrected¹⁷³ by χ^(f) = 1 and χ^(b) = e^−γτ.

In applying the LD to the construction of the TS anchored on the NHIM, the restriction of the extremization in the LD to only the forward or backward component was found useful.¹⁷³ It avoids the distortion caused by hyperbolic trajectories whose accumulations would lead to unaccountable divergences leading to unassignable LDs. If there is a manifold separating the past or future fate of trajectories, then such trajectories evolve to different regions of phase space—e.g., the reactant and product regions in a chemical reaction. In such cases, the LD is expected to describe the manifolds as a discontinuous-like abrupt change or as through a singular contour. For more details on the so-called distinguished trajectories that motivate the formulation of the LD, see Refs. 174 and 175.

The use of the LD was suggested¹²² for the construction of the TS anchored on the NHIM and implemented on a time-varying energy surface based on the Eckart barrier system in the presence of stochastic noise. The LD successfully captures the TS trajectory,^95,122 which is known as the linear approximation of the NHIM in non-autonomous systems (also known as a random fixed point and stationary orbit in the theory of stochastic dynamical systems¹⁷⁶ as summarized in the appendix of Ref. 186). When constructed in this way, the TS trajectory is denoted as the LD-TS trajectory.

LDs have also been used to find the NHIM for barrier-less reactions on a step function-like potential energy surface dissipated by a Langevin bath.¹⁶⁰ The associated time-dependent dividing surface (TDDS) is anchored at the coordinate position of the LD-TS trajectory. It successfully describes the expected monotonic population decay. Another key feature is the avoidance for the need for a reference TS, which would be needed for developing a perturbation theory but which is not available in a barrierless reaction. Nevertheless, the agreement between the LD-TS trajectory and perturbation theory for the TS trajectory was investigated in the double well potential with three different time-varying saddle points.¹⁷⁷ It was also demonstrated that the TDDS gives monotonic population decay, which does not appear if one uses a naïve dividing surface.

The LD has been implemented to uncover the reaction geometry in additional chemical systems. For example, the rates in 1D and 2D models of driven ketene isomerization reaction have been determined. The potential energy surface is represented using the Gezelter–Miller model,¹⁷⁸ and the driving is affected through coupling between a time-dependent electric field and the molecule’s dipole.¹⁶² The singular contour of the LD gives similar behavior with the boundaries of final basin at time t₀ ± τ. Patra and Keshavamurthy¹²⁷ demonstrated the similarity of the singular contour of the LD and the reactive island theory.²⁵ Specifically, they demonstrated an efficient sampling of the reaction path by noticing the similarity of the LD value for all the points in a given reactive island. Naik and Wiggins¹⁷⁹ further demonstrated the use of the LD to uncover the reactive islands in a model for conformational isomerization in a Langevin bath.

Fractal structures in the NHIM have been observed¹⁸⁰ via the singular contour of the LD. To do so, they introduced a trajectory-dependent integration time

where $q̄$ is an effective measure of TS region. Although this formulation is not strictly Lagrangian invariant, thanks to the choice of τ[q(t)], they can single-out the LD-TS trajectory. Feldmaier et al.¹⁸¹ examined the TDDS on the time-varying 2D potential energy surface using this formulation, and no recrossing was observed among the larger number—viz., 160 000— of samples that were prepared near the vicinity of the saddle point. Geometric structures in driven two-saddle systems have also been revealed using the LD.¹⁸² 

The success of the LD in enabling the construction of the NHIM and related structures showed that nonperturbative strategies could be employed for this purpose. In turn, this has led to a growing number of nonperturbative methods, including the binary contraction method,¹⁸³ the ensemble method,^184,185 local manifold analysis,^163,185 and methods based on Floquet analysis.^91,185

The asymptotic trajectory indicator (ATI)¹⁶⁴ is yet another alternative to the LD for identifying geometric structures for a reaction. It is aimed at identifying the reactivity boundaries proposed by Wall et al.⁵⁷ (and further elaborated by Wright et al.^51–54). It is based on the observation that the first passage time grows for trajectories from a given initial condition—that is within a finite area, that contains an asymptotic manifold (such as a NHIM)— as the initial condition moves closer to the stable manifold in forward time propagation. Consequently, the first passage time has a singular contour over the initial conditions at the stable or unstable manifold. It can therefore be used as an ATI to identify the location of manifold in an algorithm using neighbor bisection and continuation similar to the multidimensional bisection algorithm.¹⁸⁶ 

Many important chemical reactions occur in solution. In this section, we discuss the extension of the theory of phase space geometries for gas phase reactions into solvated reactions. The dynamical treatment of the barrier crossing problem in solution dates back to the theory of Kramers¹⁸⁷ who considered the problem through the Brownian motion of the reacting particle on a one-dimensional energy landscape with a maximum point separating the reactant and the product regions. The diffusion on the energy landscape was treated by the Fokker–Planck equation, which describes the time evolution of the probability distribution in a thermal environment. The same situation can be described in the form of an equation of motion for an ensemble of trajectories.

The Langevin equation describes the motion for each representative particle. It is expressed in the following form:¹⁸⁸ 

where q is the position coordinate for the trajectory along the reaction, V is the potential of mean force, γ is the friction constant, and ξ(t) is a random force. The random force, describing kicks from the thermal environment, is often presumed to be white noise satisfying the following statistical properties:

where k_B and T are the Boltzmann constant and the temperature, respectively. The angle brackets denote an ensemble average, and δ is the Dirac delta function.

Later, Kramers’s solution was extended by Grote and Hynes¹⁸⁹ to the situation where the dynamics is described by the generalized Langevin equation (GLE) in the following form:

where −dV/dq is the mean force, which is the equilibrium ensemble average of the force to the variable q. In this case, the friction term is given as an integration over the time from zero to the current time t with the friction kernel γ(τ). The random force in this case satisfies the following statistical properties:

The GLE has been extensively studied in the field of nonequilibrium statistical mechanics. A remarkable fact is that there is a general proof using the projection operator formulation^188,190,191 that any Hamiltonian system in thermal equilibrium can be projected onto its subsystem obeying the GLE-type equation of motion. In fact, if the “subsystem” here is taken as a set of position coordinates and their respective velocities $(q,q̇)$⁠, the most general form of the GLE-type equation of motion is written as follows:

where f(q) is the mean force at q and K in the friction term is generally a nonlinear function of the positions and the velocities in the past. The friction term is often approximated as linear in the velocities,

and the statistical property of the random force is

This form was shown¹⁹² to be exact when the bath is a collection of harmonic oscillators and the coupling to the system is bilinear. A useful consequence of this equivalent expression is the seminar result by Pollak that the Grote–Hynes rate expression cited above can be obtained using transition state theory on a rotated planar DS in the infinite dimensional space of the reactive mode and its orthogonal harmonic oscillators.^193–195 Note that the one-dimensional version of Eq. (24) is Eq. (20) and that setting γ(τ) ∝ δ(τ) results in the Langevin equation (17).

Let us start with considering the dynamics in the saddle region in the framework of the multi-dimensional Langevin equation,

where $q=(q1,q2,…,qn)T$⁠, $f(q)=(f1(q),f2(q),…,fn(q))T$⁠, and $ξ=(ξ1,ξ2,…,ξn)T$ are the column vectors of the subsystem position coordinates, the mean force, and the random force, respectively, and Γ is an n × n real matrix (γ_ij). Without loss of generality, we assume that the origin is a stationary point [f(0) = 0] corresponding to the “barrier top” on the potential. As in Sec. II, we first consider a linear approximation of the force at the origin,

Here, A is an n × n real matrix (a_ij), representing the lowest-order part of the mean force: $f(q)=Aq+Oq2$⁠. By introducing the phase space vector $(q,q̇)T$⁠, this can be written in the following form:

The matrix in the first term can be diagonalized by introducing its eigenvectors. Here, we assumed that it has one positive real eigenvalue λ₁ > 0 corresponding to the motion sliding down the barrier and the others describe stable motions $Reλj≤0for j=2,3,...$⁠. Let e_j be the eigenvectors corresponding to the eigenvalues λ_j and ${ej⋆}$ be the row vectors that are dual to {e_j},

where δ_ij is Kronecker’s delta. Phase space coordinates u_j can be introduced through the construction

or, equivalently,

Multiplying $ej⋆$ from the left to Eq. (29), we thus obtain

with

Compared to the normal mode treatment in Sec. II, the equation of motion [Eq. (34)] is more coupled in the sense that the dynamics depends on the thermally fluctuating random force. Given any initial value of the phase space coordinate u or $(q,q̇)$⁠, its future trajectory is not determined uniquely but depends on the specific instance of ξ(t). The phase space geometry as discussed in Sec. II cannot be drawn if we restrict to the phase space coordinate u that has one-to-one correspondence with $(q,q̇)$⁠. To discuss the phase space geometry for Eq. (34), Bartsch et al.^95–99 used the idea of a time-dependent shift¹⁹⁶ to define a new set of coordinates (x_j),

where the notation S is adopted from Refs. 84, 101–105, 113, and 122: for any complex value μ and a function f of time t,

Note that with this definition the following properties are satisfied:

Then, the equation of motion for x is given by

As we have assumed that the first mode is the only unstable mode (λ₁ > 0), we can construct the stable manifold that divides the phase space into two parts with respect to the future destination of the trajectories in the same fashion as explained in Sec. II. Note that the definition of the eigenvectors e_j in Eq. (30) has the ambiguity of scalar multiplication: If e_j is an eigenvector with eigenvalue λ_j, then for any arbitrary constant α, the vector αe_j is also an eigenvector with the same eigenvalue λ_j. We exploit this ambiguity to let e₁ have positive projection on the q₁ axis. Then, having x(t) → +∞ (−∞) as t → +∞ means q₁(t) → +∞ (−∞, respectively). The phase space objects are therefore as follows:

The resulting rates have also been reported in the literature.^197,198

Let us here consider a simple, one-dimensional system with an inverted parabolic potential, represented by

where $ωb2$ is the curvature of the potential at the maximum point (q = 0). The eigenvalues and eigenvectors can be calculated analytically as

The eigenvalue λ₊ and the eigenvector e₊ correspond to the unstable mode. Suppose a transition state ensemble at q = 0 with the Maxwell-Boltzmann distribution of the velocity $q̇$⁠,

The flux F is given by

where $Preaction|q,q̇$ is the probability that a trajectory with initial condition $(q,q̇)$ ends in the products’ region in the future. The original TS theory assumes $Preaction|0,q̇=1$ for $q̇>0$ and $Preaction|0,q̇=0$ for $q̇<0$⁠. Therefore,

However, we have seen that under the influence of friction and random force, the fate of a trajectory is determined by the sign of x₁ introduced above. Combining Eq. (36) and

we obtain

The reaction probability is

From Eqs. (36) and (38), with Eq. (19), we can calculate the second moment of $u1‡$⁠,

If we assume the normal distribution for the random force ξ, the probability distribution of $u1‡$ is also normal. Therefore, the probability given by Eq. (52) becomes

where the integration variable is transformed to $v=λ+−1u1‡$ to simplify the expression. Integration by parts then yields

Therefore,

This is the Kramers transmission coefficient, which was derived in Ref. 187 by the use of Fokker–Planck equation. The derivation given here demonstrates that the concept of the underlying phase space geometry provides a consistent picture.

Next, we consider the inclusion of the higher order terms in the mean force in Eq. (26). We review the treatment developed in Refs. 101–105, which is a nonlinear extension of Bartsch et al.’s work^95–99 on the linear case. By performing the power series expansion, the nonlinear equation of motion is

with the expansion coefficients α_im for each power m = (m₁, m₂, …, m_n) and monomials $qm=∏iqimi$⁠. The nonlinear part contains terms of quadratic and higher order, that is, |m| = m₁ + m₂ + ⋯ + m_n ≥ 2.

Substituting Eq. (36), the equation of motion becomes

where the coefficients f_jm(t) are obtained by substituting the transformation q ↦ x into the nonlinear part $∑|m|≥2αj,mqm$ of Eq. (60).

We employ non-Hamiltonian normal form (NF) theory¹⁹⁹ to introduce a nonlinear transformation from x to y so that the equation of motion for y₁ contains no coupling with the other coordinates up to the order of perturbation. By the transformation

the equation of motion (62) is cast into

where the two unknown functions c_j and w_j are to be determined so that the resulting equations of motion for the new coordinate y [Eq. (64)] are in as simple a form as possible.

Expanding the two functions in power series,

and performing some involved algebra,¹⁰¹ we arrive at an equation of the following form:

where $⟨λ,m⟩=def∑k=12nλkmk$ and g_jm(t) are the coefficients of the polynomial expansion of

In the process of determining c_jm(t) and w_jm(t) perturbationally order by order, g_jm(t) is known from the results of the lower orders.¹⁰¹ There thus exist two unknown quantities c_jm(t) and w_jm(t). To simplify the equation of motion for the new variable y₁, one eliminates c_1m(t) order by order by setting

using the S-symbol as in Eq. (38).

One may want to eliminate all the coupling terms between y₁ and the other coordinates so that the motion of y₁ becomes totally independent of the other DoF. Alternatively,^103,106 one can eliminate only terms with m₁ = 0 so that the equation of motion is in the following form:

where $m1′=m1−1≥0$ since the terms with m₁ = 0 have been eliminated by the normal form transformation. With the equation of motion (70), the set $Ws=(q,q̇,t)|y1=0$ is still an invariant set (i.e., dy₁/dt = 0) of codimension one. That the set $Ws$ is invariant means that $Ws$ divides the whole phase space into two parts and no trajectories can penetrate from one region to the other. A succession of time slices of the manifold along the TS-trajectory, shown in Fig. 4, illustrates the changing phase-space geometry in the reaction region and how the moving TS can be used to identify the crossing of a reactive trajectory from the reaction to product regions. As long as the linear terms (λ_jy_j) is dominant in magnitude, this division of the phase space corresponds to the classification of the trajectories with respect to their future destinations. The transformed coordinate y₁ obeying the simplified equation of motion can thus be regarded as a new reaction coordinate, while all the effects of nonlinear couplings, friction, and random forces are incorporated in the transformation from the original coordinates $(q,q̇)$ to y.

The form of Eq. (70) has been called a minimal normal form^83,103,106 in the sense that the elimination of the coupling terms has been performed to the least extent for the purpose of having $(q,q̇,t)|y1=0$ as an invariant set. It was shown to have better convergence property than other types of normal form that try to eliminate more coupling terms. See also Ref. 83 concerning the minimal normal form for Hamiltonian systems.

The discussion of the phase space structure up to here is based on the Langevin equation [Eq. (17)]. This equation can be obtained by inserting the Dirac delta function into the friction kernel in the GLE [Eq. (20)]. Physically this means that the time scale of the motion in the environment is much shorter than the subsystem. When this limit cannot be taken, one must resort to the GLE with the memory term. It may seem that the treatment of the GLE is drastically different than the LE because of the memory term involving integration as well as the differentiation. In mathematical terms, the GLE is not a usual differential equation as it takes the form of an integro-differential equation. Based on a long line of studies,^{96,105,106,196,200–207} a formulation that converts the GLE into multi-dimensional and memoryless equations of motion was recently provided.^208,209 Briefly, the GLE can be converted into the following multi-dimensional equations of motion with newly introduced dynamical variables X_m:

where α_m and c_m are coefficients of multi-exponential expression of the friction kernel,

The coefficients b_m are determined by (α_m, c_m) through the procedure detailed in Ref. 208. The new random force η(t) is a white noise,

Now that Eq. (72) is an ordinary stochastic differential equation with white noise, it can be treated with a similar manner as for the multi-dimensional Langevin equation (26). Geometric structures determining the fate of trajectories can thus be readily extracted in cases not accessible from the GLE directly.

For example, in the energy diffusion limited regime, the DS can be associated with a fixed energy—rather than position—in a transition state rate calculation.^194,195,210

The GLE and the statistical properties of the random force attached to it can be generally derived from the projection operator formulation. Usually this derivation assumes that the environment is in a thermal equilibrium state. The stationarity of the environment is reflected in the form of the friction kernel that depends only on the time difference t − t′, meaning that the way of the environment’s response to the subsystem does not depend on the time of the action from the subsystem. The equilibrium nature of the GLE is also reflected in the very fact that the statistical property Eq. (20) of the random force refers to the canonical ensemble average. On the other hand, many chemical reactions occur between non-equilibrium states. In light-induced chemical reactions, for example, the photons excite the system into a specific state (or a set of states) creating an initial distribution that are totally different from the canonical ensemble. One can even irradiate an oscillatory laser field during the whole course of the chemical reaction to guide the system into desired products. To treat the dynamics of such nonequilibrium systems, the irreversible generalized Langevin equation (iGLE) and other nonstationary versions of the GLE have been introduced in the literature.^211–220 Motivated by this work and taking advantage of the generality of the projection operator formulation, a general theory²²¹ was constructed allowing for the nonstationarity to be addressed through noise beyond the multiplicative form addressed earlier.²¹⁵ After making the linear approximation to the friction term, the theory gives the following form as the extension of Eq. (20) to the nonstationary case:

where the mean force $f(q,q̇,t)$ at time t and $q̇$ is the time-dependent ensemble average of the forces acting on the variable q. The form of Eq. (75) is called the iGLE.^214–221 

The key difference from Eq. (20) is the explicit dependence of the mean force on time t and the velocity $q̇$ and the dependence of the friction kernel γ(t, t′) on both the “initial” (t′) and the “final” (t) times. This results from the fact that the property of the bath is changing with time. The statistical property of the random force is given by²²¹