The complexity of gas and condensed phase chemical reactions has generally been uncovered either approximately through transition state theories or exactly through (analytic or computational) integration of trajectories. These approaches can be improved by recognizing that the dynamics and associated geometric structures exist in phase space, ensuring that the propagator is symplectic as in velocity-Verlet integrators and by extending the space of dividing surfaces to optimize the rate variationally, respectively. The dividing surface can be analytically or variationally optimized in phase space, not just over configuration space, to obtain more accurate rates. Thus, a phase space perspective is of primary importance in creating a deeper understanding of the geometric structure of chemical reactions. A key contribution from dynamical systems theory is the generalization of the transition state (TS) in terms of the normally hyperbolic invariant manifold (NHIM) whose geometric phase-space structure persists under perturbation. The NHIM can be regarded as an anchor of a dividing surface in phase space and it gives rise to an exact non-recrossing TS theory rate in reactions that are dominated by a single bottleneck. Here, we review recent advances of phase space geometrical structures of particular relevance to chemical reactions in the condensed phase. We also provide conjectures on the promise of these techniques toward the design and control of chemical reactions.

## I. INTRODUCTION

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.^{1} The theory of reaction rates dates back to van 't Hoff^{2} and Arrhenius^{3} 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^{34} 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^{9} and Keck^{15} 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^{20} 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,^{36} 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.^{37} 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.*^{57}), 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,^{36} 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^{114} (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.

## II. PHASE SPACE GEOMETRY IN GAS PHASE 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\u2032,q1\u2032)$ and $(p2\u2032,q2\u2032)$ are reactive and nonreactive normal form coordinates and conjugate momenta, respectively.^{130} 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\u2032,q1\u2032)$ and $(p2\u2032,q2\u2032)$ 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\u2032,q2\u2032)$ 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 2*n*-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.

### A. Perturbative analysis

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*

_{1},

where

Here, *q*_{i} and *p*_{i} are the *i*th normal coordinate and its conjugate momentum, respectively; *ω*_{i} and *C*_{jkl}, *C*_{jklm}, … are, respectively, the frequency of the *i*th 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*_{1}, *ω*_{1} ≡ − |*ω*_{1}|i (because of the negative curvature $\omega 12<0$ along *q*_{1}) and those of the other modes are pure-imaginary and real, respectively.

When the Hamiltonian *H* is the nonperturbed Hamiltonian, i.e., *H* = *H*_{0},

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\u225412\pi \u222epidqi(=12\omega i(pi2+\omega i2qi2))$ and Θ_{i} = tan^{−1}(*ω*_{i}*q*_{i}/*p*_{i}). The action and angle variables associated with the reactive DoF is formulated by $J1\u2254i2|\omega 1|(pi2\u2212|\omega 1|2qi2)$ and $\Theta 1=\u2212i\u2061tanh\u22121|\omega 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^{28} for the estimation of the tunneling reaction: $J1\u2254\u222bab2(V1(q1)\u2212E1)dq1$, where *V*_{1} and *E*_{1} are the potential energy and mode energy associated with *q*_{1}, respectively, and *a*, *b* are points satisfying *V*_{1}(*q*_{1}) = *E*_{1}.

CPT canonically transforms a coordinate system (** p**,

**) into a new coordinate system (**

*q***′,**

*p***′) 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**

*q**H*(

**,**

*p***) into a new Hamiltonian $H\u2032=H\u2032J\u2032+O(\u03f5k)$, where**

*q***′ is the new action variables associated with (**

*J***′,**

*p***′) so that $dJi\u2032/dt=O\u03f5k$ (i.e., $Ji\u2032$ is approximately time invariant), the equations of motion in the old phase space coordinates are transformed**

*q*^{64,70}into

Thus,

where the transformed new frequency $\omega i\u2032(J\u2032)$ is given by time derivative of the new angle variable $\Theta i\u2032$ conjugate to the new action variable $Ji\u2032$: $\omega i\u2032(J\u2032)=d\Theta i\u2032/dt=\u2202H\u2032(J\u2032)\u2202Ji\u2032=const+O\u03f5k$. The transformed variables $(pi\u2032,qi\u2032)$ are approximately independent of each other as long as $\omega i\u2032$ is regarded as constant in time (cf. $d2qi/dt2=\u2212\omega i2qi$ for the normal mode Hamiltonian system where each mode is decoupled) as illustrated in Fig. 3. One can derive analytical expressions of (** p**′,

**′) with respect to (**

*q***,**

*p***) and vice versa (see the full derivation in Ref. 70).**

*q*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^{32} 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\u2032(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\u2032(p,q)$ decoupled from the rest under the existence of chaos in a nonreactive space^{115} 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^{32} 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\u2032$. 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\u03f5k$ in the vicinity of saddle. The analytical expression of the reaction coordinate $q1\u2032(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**,

**) at some transit time.**

*q*### B. Analytical description of phase space geometry by canonical normal form theory

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.^{44} 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.^{83}

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^{154} 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\u03f5k$. 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.

### C. Nonperturbative analysis

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.^{166} The LDs first introduced in Ref. 167 heuristically support the conjecture that stable/unstable manifolds can be described by the “*abrupt change*”^{168} of $M$-function (e.g., arc length of a trajectory over *τ* time evolution). Mendoza and Mancho^{167} 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.^{168} Despite the controversy around the conjecture due to alleged counterexamples^{170,171} and their possible resolution^{169} and the fact that the LD is not objective^{172} (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^{110} 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*_{0} is the initial time, ** q** ($q\u0307$) is the coordinate (velocity) vector of the system,

*q*_{0}($q\u03070$) is that of the initial condition at time

*t*

_{0}, and ‖

**‖**

*x*_{p}is the

*L*

^{p}norm of vector

**defined by $\Vert x\Vert p\u2254\u2211ixip1/p$ (**

*x**L*

^{2}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

^{173}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.^{173} 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^{122} 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^{176} 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.^{160} 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.^{177} 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,^{178} and the driving is affected through coupling between a time-dependent electric field and the molecule’s dipole.^{162} The singular contour of the LD gives similar behavior with the boundaries of final basin at time *t*_{0} ± *τ*. Patra and Keshavamurthy^{127} demonstrated the similarity of the singular contour of the LD and the reactive island theory.^{25} 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^{179} 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^{180} via the singular contour of the LD. To do so, they introduced a trajectory-dependent integration time

where $q\u0304$ 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.*^{181} 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.^{182}

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,^{183} the ensemble method,^{184,185} local manifold analysis,^{163,185} and methods based on Floquet analysis.^{91,185}

The asymptotic trajectory indicator (ATI)^{164} 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.*^{57} (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.^{186}

## III. PHASE SPACE GEOMETRY IN SOLVATED REACTIONS

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^{187} 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:^{188}

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^{189} to the situation where the dynamics is described by the generalized Langevin equation (GLE) in the following form:

where −d*V*/d*q* 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\u0307)$, the most general form of the GLE-type equation of motion is written as follows:

where ** f**(

**) is the mean force at**

*q***and**

*q***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,**

*K*and the statistical property of the random force is

This form was shown^{192} 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,\u2026,qn)T$, $f(q)=(f1(q),f2(q),\u2026,fn(q))T$, and $\xi =(\xi 1,\xi 2,\u2026,\xi 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\u0307)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 *λ*_{1} > 0 corresponding to the motion sliding down the barrier and the others describe stable motions $Re\lambda j\u22640for\u2009j=2,3,...$. Let *e*_{j} be the eigenvectors corresponding to the eigenvalues *λ*_{j} and ${ej\u22c6}$ 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\u22c6$ 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\u0307)$, 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

**that has one-to-one correspondence with $(q,q\u0307)$. To discuss the phase space geometry for Eq. (34), Bartsch**

*u**et al.*

^{95–99}used the idea of a time-dependent shift

^{196}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 (*λ*_{1} > 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*_{1} have positive projection on the *q*_{1} axis. Then, having *x*(*t*) → +∞ (−∞) as *t* → +∞ means *q*_{1}(*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 $\omega 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\u0307$,

The flux *F* is given by

where $Preaction|q,q\u0307$ is the probability that a trajectory with initial condition $(q,q\u0307)$ ends in the products’ region in the future. The original TS theory assumes $Preaction|0,q\u0307=1$ for $q\u0307>0$ and $Preaction|0,q\u0307=0$ for $q\u0307<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*_{1} introduced above. Combining Eq. (36) and

we obtain

The reaction probability is

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

where the integration variable is transformed to $v=\lambda +\u22121u1\u2021$ 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.

### A. Perturbation theories

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*

_{1},

*m*

_{2}, …,

*m*

_{n}) and monomials $qm=\u220fiqimi$. The nonlinear part contains terms of quadratic and higher order, that is, |

**| =**

*m**m*

_{1}+

*m*

_{2}+ ⋯ +

*m*

_{n}≥ 2.

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

where the coefficients *f*_{jm}(*t*) are obtained by substituting the transformation ** q** ↦

**into the nonlinear part $\u2211|m|\u22652\alpha j,mqm$ of Eq. (60).**

*x*We employ non-Hamiltonian normal form (NF) theory^{199} to introduce a nonlinear transformation from ** x** to

**so that the equation of motion for**

*y**y*

_{1}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,^{101} we arrive at an equation of the following form:

where $\u27e8\lambda ,m\u27e9=def\u2211k=12n\lambda 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.^{101} There thus exist two unknown quantities *c*_{jm}(*t*) and *w*_{jm}(*t*). To simplify the equation of motion for the new variable *y*_{1}, 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*_{1} and the other coordinates so that the motion of *y*_{1} becomes totally independent of the other DoF. Alternatively,^{103,106} one can eliminate only terms with *m*_{1} = 0 so that the equation of motion is in the following form:

where $m1\u2032=m1\u22121\u22650$ since the terms with *m*_{1} = 0 have been eliminated by the normal form transformation. With the equation of motion (70), the set $Ws=(q,q\u0307,t)|y1=0$ is still an invariant set (i.e., d*y*_{1}/d*t* = 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 (*λ*_{j}*y*_{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*_{1} 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\u0307)$ 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\u0307,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}

### B. iGLE dynamics

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^{221} was constructed allowing for the nonstationarity to be addressed through noise beyond the multiplicative form addressed earlier.^{215} 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\u0307,t)$ at time *t* and $q\u0307$ 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\u0307$ 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^{221}

where $\u27e8q\u03072\u27e9t\u2032$ is the ensemble average of squared velocity at time *t*′. The average in this case is taken with respect to the given initial distribution that is not necessarily an equilibrium one.

While Eqs. (75) and (76) are proved^{221} to hold generally, the previous work^{214–220,222–224} proposed concrete forms for the nonstationary friction kernel *γ*(*t*, *t*′). These formulas provide tractable models for clear insights on the connection to and the difference from the corresponding equilibrium system. They suggested the use of a time-dependent scaling of the random force from the equilibrium one,

where the superscript “eq” denotes the (local) equilibrium state and *g*(*t*) is the scaling factor. Combining Eqs. (76) and (77) and considering the case of quasi-equilibrium^{214} between the system and the bath, $\u27e8q\u03072\u27e9t\u2032\u2248kBTbath$, where *T*_{bath} is the bath temperature, one sees that the friction kernel is also given by scaling,

A further generalization of the friction kernel was subsequently introduced^{220} allowing for multiple heat baths and time dilatation in the arguments of the equilibrium friction kernels,

where *k* labels each heat bath and *T*_{k} is the time-dependent “temperature” of the *k*th bath. The function *τ*_{k}(*t*) is a “scaled time” of the *k*th bath. This scaling of the time was originally set proportional to the temperature *T*_{k}(*t*) of the bath. This was based on the assumption that the bath evolves adiabatically with time, in which case the bath temperature was shown^{220} to be proportional to the frequency of the bath that was prescribed to vary with time. The iGLE with the friction kernel given by Eq. (79) is called the T-iGLE. In a recent discussion,^{225} the scaling was given independently from the temperature by introducing a new scaling factor *θ* as

Determining the value of the scaling factor *θ*(*t*′) by numerical fitting allows the function of Eq. (79) to adapt to a very wide range of functional forms. Reference 225 showed that a model nonstationary friction kernel, which was designed^{221} intentionally far from equilibrium, can be fitted satisfactorily to Eq. (79) by appropriately designing the scaling factor *θ*(*t*′). It is thus expected that the form of Eq. (79) is practically general and covers a wide range of function *g*(*t*, *t*′), including a system that strongly deviates from the equilibrium state.

The stationary and driven cases of the GLE have been used to characterize the rate and reaction geometry in the LiCN isomerization.^{111,163} The iGLE was used^{226} successfully to characterize the heat flows through a nanorod sliding across a surface as observed^{227} in more detailed atomistic simulations. It also provides a framework for characterizing heat flows between heat baths so as to define nonequilibrium temperatures.^{220}

The non-stationary version of the GLE has been applied to additional relevant chemical systems such as in a recent analysis of the ion dissociation of a sodium chloride ion pair (NaCl) in water.^{224} They observed damped oscillations originating from fluctuations in the water hydrogen shell. They started from the projected equation of motion of a physical variable *A*,

## IV. CONCLUSIONS AND FUTURE OUTLOOK

This paper has attempted to provide some insights in the recent development of reaction rate theory in gas and condensed phase. A modern viewpoint emerges from the analysis of the geometry and dynamics of a reaction in phase space. It originates from the spirit of Marcelin^{4,5} and Wigner,^{7,9} followed by Wall;^{57} Wright;^{51–54} and Pollak, Pechukas, and Child^{19–23} and revisited by Leon and co-workers^{25–27} in the 1980s. These developments also rely on recent advances in dynamical systems theory, which, for example, provide for invariant structures—such as the NHIM in the phase space of a many-DoF Hamiltonian—which persists even in non-Hamiltonian stochastic systems. Perhaps one of the most important points revealed by the dynamic perspectives is that one should not just elucidate the reaction rate efficiently and approximately but also unveil the origin and predictability of the fate of reactions. That is, TST would provide knowledge of the final fate of the system, even under the existence of chaos in many-DoF systems for any given initial condition in the phase space. What we have seen thus far is that this goal is possible if the reaction dynamics is dominated by a single region even if it is not addressable through perturbation theories.

We thus suggest several new directions and opportunities for resolving chemical reactions using the following space geometrical perspectives:

*Quantum manifestation of phase space geometry of reactions*: Coherent control of chemical reactions has long been a goal of physical chemistry.^{150,228–231}The recognition that the mechanism of chemical reactions is dictated by the underlying classical stable and unstable invariant manifolds provides a “lever” for such control. In particular, one could imagine a role for the light field as a “driver” for the time-dependent modulation of the invariant manifolds so as to catch and release the system from reactant to the product.^{113,114,232}However, clearly not all chemical reactions are effectively classical with respect to the atomic and molecular motion, and one should ask whether these invariant manifolds have a structure that persists in quantum regimes, which can, in turn, be used to control them. To this end, we have the success of semiclassical transition state theory^{31,32}and its recent applications to characterize gas-phase reactions.^{233–235}The underlying good quantum numbers that permit the use of this theory reflects the underlying manifolds. As such, they should be accessible for constructing quantum transition state theories that are complementary to those in recent developments.^{236–238}*Manifestation of phase space geometry of reactions in dissipative systems at equilibrium*: One might ask whether the phase space can be divided into reactive and nonreactive trajectories even for systems experiencing*stochastic*forces. Indeed, perhaps the bigger surprise is that for several cases, the answer is yes, and the better question is why. The reason partially resides in Eqs. (36) and (38): when $Re{\lambda j}>0$, one should take into account all future outcome of random force $\xi \u0303j$ into*x*_{j}at the current time*t*. In an arbitrary system in nature, this future set of random forces is unknown, and hence, such a violation of causality is not possible. In the representation provided by the Langevin equation, we can use such a formally available force sequence to effectively turn the problem into a deterministic equation and thereby use the fully specified reaction geometry to identify reactants and products. However, this is more than just mathematics. What we learn from the fluctuation–dissipation theorem is the statistical properties of random force, and it is satisfied by the forces from nature, not just the Langevin representation. If the analytical formula of the actual reaction coordinate exists under the existence of chaos and stochastic forces, it provides us with the*analytical*expression of an exact reaction probability of how many trajectories can end up with the product on average under infinitely many realization of the random force.^{97,101,102}How would the existence of such analytical expression be related to efficiency of biological functions in complex environment? And how do we begin to construct the time-dependent structure of the reaction geometry—viz, the TS trajectory—when we do not know the force sequence from such a complex environment?*Manifestation of phase space geometry of reactions in dissipative systems at nonequilibrium*: The iGLE provides a general platform for dissipative systems in which one does not need to assume canonical or stationary distributions around the reaction center. This is the more typical scenario for reactions that are driven,^{185,239,240}for example, by a laser excitation of chemical bonds near the reaction center (see, e.g., Ref. 241). How can one obtain the corresponding reaction rate in an arbitrary nonequilibrium and nonstationary environment?*Rare event sampling vs phase space geometry*: From the outset of the work on the development of the geometric perspective on transition states for dissipative systems, it was recognized^{95,99}that the geometry of the dividing surface is complementary to the ensemble in transition path sampling.^{124,125}Recently, this close relationship was further addressed by connecting the results of transition path sampling and the underlying reactive tubes.^{127,242}An open question remains as to how can the concept and tools of of the phase space reaction geometry be utilized to enhance rare event sampling.*Data scientific approaches to extract phase space geometry of reactions*: Machine learning and more general data science approaches are currently being implemented widely to address the challenges of complexity and multi-dimensionality that also plagues attempt to obtain reaction rates and pathways. One increasingly common approach is the use of machine learning to accelerate the integration of trajectories. While this is useful, in principle, it does not necessarily take advantage of the reaction geometry or enable its determination. Nevertheless, it is possible to use machine learning to extrapolate the NHIM from a dataset of high accuracy points determined by methods such as those discussed here^{243,244}or even to discover these structures from trajectories.^{245,246}The question remains as to how to better leverage machine learning tools to determine the reaction geometry of complex and high-dimensional reactions.

## ACKNOWLEDGMENTS

We thank Professor Shinnosuke Kawai for computing the second order normal form coordinates and velocities in Fig. 4 and various valuable comments and suggestions on this manuscript. We deeply acknowledge all our collaborators and colleagues working on this interdisciplinary subject between reaction dynamics and their phase space geometry and all pioneers studying transition state theories along dynamic (phase space) perspectives. This work was partially supported by the JSPS, the Cooperative Research Program of “Network Joint Research Center for Materials and Devices,” Grant-in-Aid for challenging Exploratory Research, Grant-in-Aid for Scientific Research (B), and the National Science Foundation (NSF) through Grant No. CHE 1700749.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

All the authors wrote the manuscript equally. Y. N. created all the figures. R.H. and T.K. are equally responsible for correspondence.

## DATA AVAILABILITY

Data sharing is not applicable to this article as no new data were created or analyzed in this study.

## REFERENCES

*Studies in Chemical Dynamics*(London, 1896).