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.

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.

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 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.

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 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.

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