Dynamics and decay rates of a time-dependent two-saddle system
Johannes Reiff, Matthias Feldmaier, Jörg Main, Rigoberto Hernandez
DDynamics and decay rates of a time-dependent two-saddle system
Johannes Reiff, Matthias Feldmaier, J¨org Main, and Rigoberto Hernandez
2, 3, ∗ Institut f¨ur Theoretische Physik I, Universit¨at Stuttgart, 70550 Stuttgart, Germany Department of Chemistry, Johns Hopkins University, Baltimore, Maryland 21218, United States Departments of Chemical & Biomolecular Engineering, and Materials Science and Engineering,Johns Hopkins University, Baltimore, Maryland 21218, United States (Dated: February 26, 2021)The framework of transition state theory (TST) provides a powerful way for analyzing the dy-namics of physical and chemical reactions. While TST has already been successfully used to obtainreaction rates for systems with a single time-dependent saddle point, multiple driven saddles haveproven challenging because of their fractal-like phase space structure. This paper presents the con-struction of an approximately recrossing-free dividing surface based on the normally hyperbolicinvariant manifold in a time-dependent two-saddle model system. Based on this, multiple meth-ods for obtaining instantaneous (time-resolved) decay rates of the underlying activated complex arepresented and their results discussed.
I. INTRODUCTION
One of the central aims in the field of chemical reac-tion dynamics is the accurate determination of reactionrates. This is not just an abstract problem of academic(or basic) concern [1–4], but also a practical problem withmany potential applications in complex reactions [5–9].The possibility of optimizing reaction rates by externaldriving could perhaps take these applications further inoffering improvements to throughput and efficiency.Multi-barrier reactions were considered early [10, 11]in the context of quantum mechanical tunneling throughbarriers at constant energy. In this
M-problem , the peri-odic orbit between the barriers gives rise to the possibilityof an infinite number of returns to the turning point fromwhich tunneling can proceed. The return times are usu-ally not commensurate with the period, either because ofcoupling to other degrees of freedom—such as from thebath—or because of variations in the potential. In suchcases, the coherence in the returns is altered, changingthe nature of the dynamics in ways that we address inthis work.Problems involving fluctuating [12, 13] or oscillat-ing [14–16] barriers have also received significant at-tention leading to, for example, the identification ofthe phenomenon of resonant activation [17, 18]. Whilethe approaches originally focused on the overdampedregime [12, 14], underdamped systems were later exam-ined [16, 19]. For example, mean first passage timeshave been employed to calculate (diffusion) rates in spa-tially periodic multi-barrier potentials. Therein, vari-ous static [20] as well as stochastically driven [21, 22]cases have been characterized primarily through numer-ical methods.In the current context, the main challenge in a multi-saddle system comes from the unpredictability of states ∗ Correspondence to: [email protected] in the intermediate basin. A reactant entering this re-gion may leave either as a reactant or product depend-ing on the exact initial conditions [23–26]. Historically,this challenge has been approached by categorizing reac-tions into two classes [24, 27, 28]:
Direct reactions ex-hibit a single transition state (TS).
Complex reactions,on the other hand, have two clearly separated TSs. Thepotential well between those barriers is assumed to besufficiently deep that it gives rise to a long-lived colli-sion complex. Trajectories passing through one TS enterthis collision complex and hence cannot be correlated totrajectories passing through the other TS. In reality,however, a reaction cannot always be uniquely classified.These concerns were addressed in a unified theory byMiller [24], and later refined by Pollak and Pechukas [29]so as to address shallow potential wells. While an impor-tant advancement, this theory still treats the saddle’s in-teractions statistically, thereby neglecting dynamical ef-fects like resonances. Moreover, they considered multi-step reactions in which the positions and heights of thebarriers are time-independent. Last, there are numerouspublications on valley-ridge inflection points, which aretypically described by a normal TS followed by a shared one [26, 30, 31].Craven and Hernandez [32] recently examined a four-saddle model of ketene isomerization influenced by atime-dependent external field. They encountered compli-cated phase space structures similar to those in systemswith closed reactant or product basins [33]. As a result,their analysis was limited to local dividing surfaces (DSs)and no reaction rates were calculated. Moreover, success-fully calculating instantaneous rates based on a globallyrecrossing-free DS attached to the normally hyperbolicinvariant manifold (NHIM) of a time-dependent multi-saddle system has—to our knowledge—not yet been re-ported.In this paper, we address the challenge of determin-ing the instantaneous TS decay rate for systems that notonly feature multiple barriers along the reaction path,but that are also time-dependent. Specifically, we investi- a r X i v : . [ phy s i c s . c h e m - ph ] F e b gate an open, time-dependent two-saddle model with onedegree of freedom (DoF) as introduced in Sec. II A. Thetheoretical framework along with the numerical methodsused throughout the paper are described in Secs. II Bthrough II D. We then discuss the system’s phase spacestructure under the influence of periodic driving at dif-ferent frequencies in Sec. III. By leveraging unstable tra-jectories bound between the saddles, we can construct analmost globally recrossing-free DS associated with the so-called geometric cross in Sec. IV. The DS is used in Sec. Vto calculate instantaneous decay rates and averaged rateconstants for the underlying activated complex using dif-ferent methods. Thus we report a detailed analysis ofthe nontrivial phase space structure of the chemical M-problem and the calculation of its associated decay rates. II. METHODS AND MATERIALSA. Time-dependent two-saddle model system
In this paper, we investigate the properties of multi-barrier systems by considering a 1-DoF model potentialfeaturing two Gaussian barriers whose saddle points arecentered at x = ±
1. Initially, both barriers are placedat the same level. As we are interested in consideringthe time-dependent case, however, we drive the barrier’sheights B ϕ ( t ) sinusoidally in opposite phases. That is, weuse the same amplitude and frequency ω for both saddles,but opposite initial phases ϕ ∈ { , π } . This leads to thepotential V ( x, t ) = B ( t )e − ( x +1) + B π ( t )e − ( x − (1a)with B ϕ ( t ) = 74 + 14 sin( ωt + ϕ ) . (1b)The oscillation frequency ω is a free parameter that canbe varied relative to the other natural time scales of thesystem at a given fixed total energy. At an arbitrarytime, one of the barriers will be larger than the other.For example, when the second barrier is larger, the po-tential takes a shape such as that shown at the top ofFig. 1(b). Throughout this paper, dimensionless unitsare used to explore the range of phenomena that canarise from varying the relative timescales of the systemand the driving. B. Dividing surfaces and the NHIM
Reactants and products on the canonical potential en-ergy surface (PES) are usually separated by a DS of-ten associated with a rank-1 saddle, such as that shownin Fig. 1(a), whose unstable direction identifies the re-action coordinate. The maximum energy configurationalong the corresponding one-dimensional minimum en-ergy path [34–36] x is called the TS [3, 28, 37, 38]. In V (a) DSR P (b) local DS global DSR I P x v x W s W s W u W u x FIG. 1. Typical structures of static potentials V ( x ) and theircorresponding phase spaces x – v x with one and two barriers inpanels (a) and (b), respectively. The potential barriers sep-arate reactant (R) and product (P) states. The two-barriercase features an additional intermediate (I) state in between.The maxima are associated with a hyperbolic fixed point (di-amonds) and a dividing surface (indicated by dashed verticallines) each. The corresponding manifolds divide the phasespace into four distinct, numbered regions in panel (a) andsix regions in panel (b). See Secs. II and III A for details. transition state theory (TST), the local decay rate is ob-tained from the flux through the DS.A DS is locally recrossing-free if no particle pierces itmore than once before leaving some pre-determined inter-action region around the saddle. In this case, the decayrates (to exit the interaction region) as determined by theDS are locally exact. We refer to a DS as being globallyexact or recrossing-free if the above is true independentof the choice of the interaction region as long as said re-gions do not overlap with the stable reactant or productregions [39]. Using this definition avoids inherent recross-ings caused by reflections in closed reactant or productbasins [33]. In this paper, we only address transitionsover barriers in series, and we do not address the paral-lel case in which a reaction could access more than onedistant barrier. The scope of the definitions of globallyexact and recrossing-free is therefore limited accordingly.Every saddle point of a d -DoF potential is associatedwith (2 d − W s and W u in phase space—cf. Fig. 1(a). Their (2 d − d − C. Revealing geometric structures
The geometric structure of the phase space can be re-vealed using the Lagrangian descriptor (LD) [32, 49–52]defined by L ( x , v , t ) = (cid:90) t + τt − τ d t (cid:107) v ( t ) (cid:107) (2)for a given initial position x , velocity v , and time t . Itmeasures the arc length of a trajectory x ( t ) in the timeinterval t − τ ≤ t ≤ t + τ . A local minimum in the LDarises when the particle covers the minimum distance inthe interval t − τ ≤ t ≤ t + τ . It consequently re-mains longer in the interaction region when integratingforwards or backwards in time, and thus provides a sig-nature for the presence of a stable or unstable manifold,respectively.The LD has the advantage that it is conceptually verysimple and that it can be applied to practically any sys-tem. This makes it suitable for a first visual inspection.As discussed in Ref. [53], however, it features a nontriv-ial internal structure. Numerically determining the exactposition of stable and unstable manifolds is therefore dif-ficult.A numerically simpler scheme is based on the con-cept of reactive (and nonreactive) regions as describedin Refs. [39, 53]. It discriminates initial conditions byfirst defining an interaction region in position space thatencompasses the relevant dynamics. Particles are thenpropagated forwards and backwards in time until theyleave said interaction region. In both directions of time,a particle can end up as either reactant (R) or product(P). This leads to four possible classifications for a giveninitial condition as shown in Fig. 1(a) for each of thefour regions: 1. nonreactive reactants R → R, 2. nonre-active products P → P, 3. reactive reactants R → P, and 4. reactive products P → R. Similar concepts have beenintroduced in, e. g., Ref. [32].Such a classification for the initial conditions has to beextended to include the consequences of a local minimumbetween the saddles of the reacting system of Eq. (1). Alow-energy particle trapped near this local minimum [cf.Fig. 1(b)] would lead to diverging computation times be-cause it may never leave the interaction region. To solvethis problem, an additional termination condition is in-troduced, whereby any particle that crosses the potentialminimum a specified number of times n mc is classified asan intermediate (I) particle. Consequently, up to ninedifferent regions in phase space can be distinguished forany given value of n mc .Using the concept of reactive regions, stable and un-stable manifolds can be revealed as borders between ad-jacent regions. Their closure’s intersections form theNHIM. The algorithm used to calculate points of theNHIM is based on the binary contraction method (BCM)introduced in Ref. [53]. It starts by defining a quadran-gle with one corner in each of the phase space regionssurrounding the manifold intersection. The quadrangleis then contracted by successively determining the regioncorresponding to an edge’s midpoint and moving the ap-propriate adjacent corner there. The BCM is thereforeeffectively composed of four intertwined classical bisec-tion algorithms. To reliably identify the initial corners,we modify the BCM slightly as detailed in Appendix A. D. Calculating decay rates
The existence of a NHIM of codimension 2 and its rolein determining the chemical reaction rate brings an ad-ditional concern. Namely, what is the degree of insta-bility of the TS as determined by the decay of trajec-tories that start in the proximity of the NHIM? In atime-dependent—e. g. driven—environment, this instan-taneous decay will be time-dependent as well. Neverthe-less, it can be assigned a single characteristic decay rateconstant when the time-dependence is periodic by takingthe average over the period [14, 54].In this paper, we implement three different meth-ods, summarized in Appendix B, for calculating decayrates: (i)
The ensemble method [14, 54] yields instanta-neous (time-resolved) rates by propagating a large num-ber of particles. It is computationally expensive butconceptually simple. (ii)
The local manifold analysis(LMA) [54, 55] accelerates the computation of instan-taneous rates by leveraging the linearized dynamics nearthe NHIM. (iii)
If only average rate constants are de-sired, the Floquet rate method [54, 56] can be used in-stead while requiring even less computational resources.In the cases resolving the dynamics of the two-barrierproblem of Eq. (1), all three generally converge withinreasonable time. However, they each involve different as-sumptions which might have led to different results, andwhich can provide complementary interpretations aboutthe underlying dynamics. As shown in the results sec-tions, all three lead to decay rates in excellent agreement.The repetition thus also serves to provide assurance inthe reported values.
III. GEOMETRIC STRUCTURE OF THETWO-SADDLE SYSTEM
The phase space structure of the model system intro-duced in Sec. II A is highly dependent on the drivingfrequency ω . In the following we will give a qualitativeoverview of the behavior the system can exhibit. A. Limiting cases
A static two-saddle system akin to Eq. (1) with ω = 0exhibits the phase space structure shown in Fig. 1(b).The saddle tops are associated with hyperbolic fixedpoints whose stable and unstable manifolds each forma cross. If the first saddle is smaller than the second one,two of its manifolds constitute a homoclinic orbit. In thiscase, the phase space is composed of six regions, namely1. nonreactive reactants R → R, 2. nonreactive productsP → P, 3. reactive reactants R → P, 4. reactive productsP → R, 5. particles that react over the first saddle butget reflected at the second R → I → R, and 6. inter-mediate particles that are trapped between the saddlesI → I.Likewise, if the driving frequency is sufficiently large( ω → ∞ ), the particle will effectively see an averagestatic potential in which it must cross two similar staticbarriers of equal height. As the energy is conserved, oncethe particle crosses the first barrier, it necessarily crossesthe second barrier. This results in a phase space structuresimilar to that for the static case of Fig. 1(b)—althoughwith heteroclinic orbits connecting the hyperbolic fixedpoints.In such (effectively) static cases, it is straightforwardto define a global recrossing-free DS. In a 1-DoF con-stant energy system, if (and only if) a particle crossesthe highest saddle, it has demonstrated to have enoughenergy to react over all saddles. Since the largest barriertherefore unambiguously determines whether a particlereacts or not, its associated local DS becomes the global(recrossing-free) DS.While this holds true for static systems, dynamicallydriven systems may exhibit much richer dynamics. Forexample, in the case of an alternating pair of dominantbarriers, as in the model of Eq. (1), the naive DS jumpsdiscontinuously from one side to the other twice per pe-riod. As a result there exist reactive trajectories thatnever cross the DS. Instead, the latter jumps over theformer leading to an inconsistent description of the reac-tion. Addressing this issue by defining particles betweenthe local DSs as reacting the moment the dominant sad-dle changes would lead to more problems, e. g., unphysi- − x − v x (a) min max L (b) × (c) × (d) × FIG. 2. Phase space structure for the time-dependent po-tential in Eq. (1) with ω = π at t = 3 / L given in Eq. (2) with τ = 16.Although the two geometric crosses seen in Fig. 1(b) are stillpresent in panel (a), they now exhibit a complicated substruc-ture involving a vast number of homoclinic and heteroclinicpoints as well as homoclinic and heteroclinic orbits. The pro-gressively zoomed cutouts shown in subpanels (b) through (d)exhibit self-recurring structures. Labels in the bottom rightcorners indicate the corresponding enlargement from the pre-vious zoom level. The potential at time t is indicated in thetop left inset of panel (a). cal Dirac delta peaks in the reaction rate. To solve thisissue, the system has to be treated as a whole. B. Intermediate driving frequency
The limiting cases discussed so far result in effectivelystatic systems. Since we are interested in novel and non-trivial behavior, however, we will now turn to interme-diate driving frequencies. These can exhibit varying de-grees of complexity as a function of the driving frequency ω . An example of such non-trivial behavior with a highlycomplex phase space is shown in Fig. 2. We use the LDdefined in Sec. II C for visualization since it is very gen-eral and requires little knowledge about the system.The geometric structure was obtained for the drivenpotential V ( x, t ) of Eq. (1) at an intermediate drivingfrequency ω = π . The general shape of the boundariesseparating reactive and nonreactive regions (cf. Fig. 1(b))is still vaguely visible. However, the precise position ofthe crossing points between the stable and unstable man-ifolds can no longer be determined. This family of cross-ing points together with the associated stable and un-stable manifolds within their vicinity appears as a crossthat has arisen from all of these geometric considerations.For simplicity, we define it as a geometric cross through-out this work. Note that this term is not meant to be aprecise mathematical structure but rather an illustrativeconcept for describing the complex phase space of thesystem under study.The fractal-like geometric crosses seen in the series ofFigs. 2(a)–(d) for a finite τ suggests a fractal structureat all scales for τ → ∞ . This structure emerges fromparticles trapped between the two saddles. For example,a reactant can enter the intermediate region over the leftsaddle, be reflected multiple times at both saddles, andfinally leave the interaction region over the right saddleas a product. The number of reflections is in this casehighly dependent on the particular initial conditions asa result of the system being chaotic. In turn, this leadsto a discontinuity in the LD, and eventually to the self-recurring patterns of a fractal structure.Figure 2 also supports the observation of geometricstructures at a given time that are thinner near the dom-inant saddle compared to the lower energy saddle. InFig. 2(a), when the right barrier is dominant, particlesinitialized near the higher saddle start with higher po-tential energy and therefore have a lower chance of be-ing reflected. As a result, fewer of these particles lingerin the interaction region and the phase space structuresare thinner. While this may be an interesting observa-tion, the geometric cross near the dominant barrier isstill highly fractal. There are no isolated, weakly frac-tal geometric crosses that could reasonably be trackednumerically. Consequently, we cannot make meaningfulstatements about a globally recrossing-free DS. C. Slow driving frequency
The previous sections have shown the range of com-plexity the model system (1) can exhibit. We now needto move from the aesthetically pleasing structures ofFig. 2 to a more rigorous identification of the globallyrecrossing-free DS. To do so, we switch to a lower drivingfrequency ω = π/
10 which is simpler to analyze but stillexhibits nontrivial behavior. Additionally, we employ theconcept of reactive regions as described in Sec. II C in-stead of the LD. The partitioning of the phase space intonine distinct regions allows to make quantitative assess-ments more easily.Application of this analysis to V ( x, t ) with ω = π/ → P on top, R → R to the left, P → P to the right,and P → R below. In the following, we refer to thisgeometric cross as the primary geometric cross . IV. THE GEOMETRIC CROSS
We will now analyze the primary geometric cross andits associated TS trajectory in more detail. Its time-dependent position can be tracked precisely over a fullperiod using the algorithm described in Sec. II C and Ap-pendix A. The result is marked as a series of black dotsin Fig. 3. The corresponding trajectory connecting thosedots is indicated as a black dashed line.
A. Global transition state trajectory
The primary geometric cross associated with the in-stantaneously higher saddle in Fig. 3 remains on the bar-rier nearly as long as the barrier remains dominant. How-ever, when the barriers’ heights approach each other, thegeometric cross quickly moves from one barrier to theother in the following way. The geometric cross beginsto rapidly accelerate towards the middle ( x = 0). Itcrosses the local potential minimum exactly when bothsaddle points are level (e. g. t = 10) and continues inthe same direction until it is located near the now highersaddle (e. g. t = 15). The reverse happens in the follow-ing half period, thereby forming a closed trajectory withthe same period T = 20 as the potential V ( x, t ) (dashedline in Fig. 3).This primary geometric cross marks the position of aparticle on an unstable periodic trajectory trapped in theinteraction region. The particle on this trajectory oscil-lates between the saddles with a period of T = 20 sothat it is always located near the higher saddle. Manyother unstable periodic trajectories associated with ge-ometric crosses—i. e. hyperbolic fixed points—in phasespace have also been found for particles in the interactionregion. Two such trajectories are shown in Fig. 4. Theyare typical of an increased degree of structure relative tothe primary TS trajectory. All such trajectories togetherform the system’s disjoint, fractal-like time-dependentNHIM. In contrast to the primary TS trajectory, how-ever, all other trajectories have periods larger than T .The only exceptions to this observation are the local TStrajectories in the vicinities of the saddle maxima, whichare also part of the global NHIM. These trajectories havethe same period T as the potential by construction.The complex dynamics that arises from the movingbarriers is also affected by the degree to which a given tra-jectory is decoupled from the barriers as it traverses thewell between them. Although the energy of the minimum − v x t = 0 2 . . − x − v x . − x − x . − x → R I → R P → R R → I I → I P → I R → P I → P P → P FIG. 3. Reactive regions for potential (1) with ω = π/
10 and n mc = 8 as a function of time t (cf. inset labels). The color (orshaded) legend at the top indicates whether a particle starts or ends in the reactant (R), intermediate (I), or product (P) state.The position of the point on the NHIM associated with the primary geometric cross (black dot) is tracked across a full period(TS trajectory, black dashed line). The insets in the upper right corner of the top row panels and the lower right corner of thebottom row panels illustrate the potential at the corresponding initial time t . − − −
20 0 20 40 60 t − x FIG. 4. Examples of typical periodic trajectories x ( t ) withperiods T = 40 (darker gray or blue) and T = 60 (lightergray or orange) for potential (1) with ω = π/
10. The primaryTS trajectory from Fig. 3 with period T = 20 (black dashedline) is shown for comparison. stays roughly constant, as indicated in Fig. 3, its positionmoves back and forth between the barriers as shown inFig. 5(a) so that it is always closer to the lower barrier.This is a manifestation of the Hammond postulate [57]but now applied to the negative potential. When the first (left) barrier dominates the dynamics, a locally non-recrossing DS associated with it identifies the trajectorieswith sufficient energy to cross it, and those continue un-abated across the well and the second barrier. As a con-sequence, a DS located at the well identifies the reactivetrajectories equally well (or badly) in this regime. Thisis indicated in Fig. 5(b) by good identification when theactivated particle continue past the second barrier or dis-agreement when misidentification of trajectories begin tobe reflected across both. When the second (right) barrierdominates the dynamics, identifying reactive trajectoriesat a DS at the distant well allows the evolving trajecto-ries to be reflected by the barrier leading to recrossings.Thus we find that the reaction dynamics in between thebarriers—e. g. at the well—does not go through a singleidentifiable doorway. In turn, this points to the need fordescribing the dynamics—even in a local sense—througha geometric picture that spans the two barriers, e. g. theglobal NHIM.Meanwhile, since the NHIM now consists of more thanone trajectory, it is not necessarily obvious which of theseis most suited for attaching a global DS. We can, how-ever, set conditions the global TS trajectory should fulfill.First, for symmetry reasons, the trajectory should havethe same period T = 20 as the potential. Second, the − x D S (a) t e rr o r [ % ] (b) local (left)local (right)discontinuouspot. min.TS traj. FIG. 5. (a) Position x DS and (b) percentage of trajectorieswith recrossings or classification errors as a function of time t for various choices of their assignment as reactive or nonreac-tive. The assignment is performed according to the crossingof a specified DS associated with either the local (solid) orglobal (thick dashed) TS trajectories, the instantaneous po-tential minimum between the saddles (dotted), or the discon-tinuous TS trajectory jumping between the local ones (dash-dotted). For every t , an ensemble of 10 particles with uni-formly distributed velocities 1 . ≤ v x ( t ) ≤ . x ( t ) = − t = 80 time units. Theattached DS is parallel to v x for all times. The right local DSis recrossing-free by construction since particles are startedsolely to the left of the saddles and is therefore not shown in(b). This specific choice for the ensemble is also the reasonfor the graph’s asymmetry. global TS trajectory should approach the higher saddle’slocal TS trajectory in the (quasi-)static limit ω → B. Comparison of dividing surfaces
The next task is to determine the degree to which theglobal TS trajectory gives rise to a recrossing-free DS.Numerically, this can be tested by attaching a DS to itas specified in the caption of Fig. 5, propagating an en-semble initialized near it, and recording the number and direction of DS crossings (or not) that transpire there-after in the propagated trajectories. For simplicity, weconsider only the most challenging cases in which theensemble’s initial energies are chosen to be between thesaddles’ minimum and maximum heights. The usual er-ror in the DS is signaled by the existence of more thanone crossing for the trajectories, and the fraction of suchrecrossings is used below as a measure for the DS’s qual-ity.For simplicity, we limit ourselves to DSs defined by x = x DS ( t ), i. e., parallel to v x . The results of this analy-sis are shown in Fig. 5(b). As can be seen, the global DSassociated with the time-dependent geometric cross fea-tures error rates that are significantly reduced comparedto the local DS fixed at the left barrier. One possibleDS could be constructed by placing it at the instanta-neous potential minimum—shown as the dotted curve inFig. 5(a). It would be expected to be ineffective giventhat rates are usually determined by rate-limiting barri-ers, not valleys, in between reactants and products. In-deed, the recrossing errors found for this DS, shown inFig. 5(b), were high and even worse than those from theuse of the DS fixed at the left barrier. But the highesterror rate comes from the naive attempt to treat the DSassociated with the instantaneously higher saddle pointas the global one. In this case, errors can arise fromevents beyond the recrossing of the DS. That is, therenow exists the possibility that the discontinuous instan-taneous DS can jump over the trajectory. It is the com-bination of recrossing and classification errors that leadsto the jagged and large deviations in the % error seen forthe naive discontinuous DS. As discussed in Appendix C,this can lead to misclassification of the reactivity of thetrajectory.Finally, we consider a local DS fixed at the right bar-rier. This choice would lead to no recrossing or classifica-tion errors for trajectories moving in the forward direc-tion (from reactants to products). However, it would failbadly for trajectories moving in the backward directionby symmetry with our finding for the forward trajecto-ries crossing the DS at the left barrier. Thus, the useof the global DS associated with the TS trajectory bestcaptures the time-dependent geometry of the reaction.Although this TS-trajectory DS is much better thanany alternatives considered so far, it still exhibits anamount of recrossings that cannot be explained bynumerical imprecision alone. Instead, recrossings arecaused by the fractal-like phase space structure of thesystem: Figure 6(a), for example, shows the phase spacestructure at t = 0 .
25. We can see a relatively large patchof nonreactive reactants (labeled R → R, light blue) tothe right of the DS. Particles in this patch leave the in-teraction region to the reactant (left) side forwards andbackwards in time. Consequently, they need to cross theDS at least twice, which counts as a recrossing. An anal-ogous argument can be applied to P → P regions to theleft of the DS. The fractal-like phase space structure thusleads to numerous problematic patches of vastly different − . . x − v x (a) − .
001 0 x − x DS ∆ x ∆ v u x W u W s (b) × × FIG. 6. Reactive regions analogous to Fig. 3 at t = 0 .
25 withthe addition of the DS (solid black line). (a) The simplestchoice of a DS parallel to v x leads inevitably to recrossingsas indicated. (b) The immediate region of the TS trajectory[shown as a dot in panel (a)], is enlarged, as indicated, toreveal the geometry of the ensemble (densely dotted) usedin the rate calculations of Sec. V. The ensemble is sampledequidistantly on a line parallel to the unstable manifold W u at distance x − x DS = − × − from the TS trajectory at x DS ( t ). The manifolds W u and W s are given by the bound-aries between reactive and nonreactive regions. The differen-tial of W u is indicated by the sides ∆ x and ∆ v u x of the slope. sizes, hence the two distinct peaks in Fig. 5(b).A totally recrossing-free DS, by contrast, would neces-sarily have to divide the phase space such that R → Rregions are always on the reactant side, and P → P re-gions always on the product side. This is not possiblewith a planar DS of any orientation due to the system’sfractal-like nature. A globally recrossing-free DS—if itexists—would have to curve as time passes because theentirety of the phase space between the saddle points hasa clockwise rotating structure (including periodic trajec-tories on the NHIM, cf. Figs. 3 and 4). The periodicityof this system would then result in a fractal, spiral-likeDS. Thus the next step in generalizing this theory wouldrequire the identification of a non-planar DS anchored atthe TS trajectory which we leave as a challenge to futurework.
V. REACTION RATE CONSTANTS
We can now calculate decay rate constants k for theactivated complex using the global TS-trajectory DS de-fined in Sec. IV. The patches leading to recrossing (cf.Sec. IV B) are unproblematic in this case because theensemble was selected close to the NHIM and therebynecessarily far from them, as can be seen in Fig. 6. Thefew recrossings that do still occur are artifacts from thenumerical error in the propagation, and are sufficientlysmall in number that their effect is smaller than the nu- t k ensemblemanifoldsFloquetlocal FIG. 7. Various instantaneous reaction rates k parameterizedby time t and associated with the global TS trajectory of thepotential defined in Eq. (1) with ω = π/
10, and driving period T = 20. The ensemble rate k e ( t ) (blue solid) is obtained bypropagation of 20 ensembles of 10 particles each as describedin Appendix B 1. The manifold rate k m ( t ) (orange dashed) isobtained using Eq. (B2). The Floquet rate constant k F (greendotted) is obtained from Eq. (B5). The mean rates, ¯ k e and¯ k m (thick horizontal lines) are averaged over the period asdefined in Eq. (3). For comparison, the Floquet rate constant k locF (red dash-dotted) of a single barrier’s local TS trajectoryis also shown. merical precision of the calculation.In the following, we consider the three different pos-sible approaches defined in Sec. II D in order to demon-strate their equivalence in multi-saddle systems. An ex-ample for the initial reactant ensemble and manifold ge-ometry can be found in Fig. 6(b).While the application of the LMA and the Floquetmethod to our model system are straightforward, ap-plying the ensemble method poses a challenge: A finiteensemble of reactants—e. g. of size 10 as implementedhere—will mostly react within a short time—e. g. 2 to 5units of time in the case shown here—compared to theperiod of driving T = 20. Resolving the whole periodwith a single ensemble would therefore require an ex-ponentially growing ensemble size. This would not benumerically feasible. Instead, multiple ensembles havebeen started at times t j incremented at equal intervals∆ t , and set to 1 in the current case. An instantaneousrate k e ( t ; t j ) can be obtained for each ensemble j . The in-stantaneous rate for the whole period k e ( t )—independentof t j —can then be recovered by concatenating the seg-ments of each k e ( t ; t j ) for t from t j to t j +1 . Comparedto the first option, this approach scales linearly with thesystem’s period instead of exponentially leading to vastlydecreased computing times and increased numerical sta- TABLE I. Values of the reaction rate constants discussed inFig. 7. The averaged ensemble and LMA rate constants ¯ k e and ¯ k m match the Floquet rate constants to within less than0 . k e . k m . k F . k locF . bility.The results are shown in Fig. 7. Since k e ( t ) and k m ( t )are not constant in time, we additionally show the aver-age ¯ k = 1 T (cid:90) T d t k ( t ) (3)over one period T of the TS trajectory.Figure 7 shows a large variation in the instantaneousreaction rate in the interval 0 . (cid:46) k ( t ) (cid:46) .
85. Twofeatures are distinctive. First, the rate k ( t ) shows mostlyflat plateaus during the time the TS trajectory is locatednear a saddle. This can be explained using Fig. 3: Al-though the phase space structure as a whole is undergo-ing significant changes during, e. g., 2 . ≤ t ≤ .
5, thelocal vicinity around the TS trajectory stays almost un-affected. Second, there are deep dips in k ( t ) while the TStrajectory moves between the saddles (as seen at around t ∈ { j | j ∈ Z } ). As can be seen in Figs. 3 and 6, thesetimes are characterized by much more shallow geometriccrosses with a low difference in the slopes of the stableand unstable manifolds. This effect is particularly appar-ent in Fig. 6(b).The same observations can also be interpreted anotherway. By comparing Fig. 7 to Fig. 5(a), we can see a clearcorrelation between the velocity of the TS trajectory andthe instantaneous rate k ( t ): the faster the TS trajectorymoves, the lower the rate drops.As can be seen in Fig. 7 and Table I, ¯ k e , ¯ k m , and k F arein excellent agreement, which illustrates the equivalenceof all three methods. The local Floquet rate constant k locF of a single saddle, on the other hand, differs significantlyfrom k F , even though both saddles are identical. Whilethe local rate constant can thus be used as an upper limitfor the overall rate constant, there is no straightforwardway to derive a global rate constant from it. Thus globalmethods employing the full TS trajectory are necessary ifaccurate rates for multi-saddle systems are desired. Allthree methods in this section satisfy this requirement,and are consequently in agreement. VI. CONCLUSION AND OUTLOOK
In this paper, we have fully characterized the reactiongeometry and determined the associated decay rates inoscillatory (or time-dependent) two-saddle systems in theframework of TST.The first set of central results of this work lies in re-vealing the phase space structure of the two-saddle modelsystem. While the structure of stable and unstable mani-folds is straightforward for (quasi-)static ( ω →
0) or veryfast oscillating ( ω → ∞ ) systems, intermediate frequen-cies lead to a fractal-like phase space. In this case, the ex-istence of a completely recrossing-free DS is questionable.For lower oscillation frequencies, however, an isolated ge-ometric cross with negligible substructure—referred to asthe primary geometric cross—emerges. This structure ispart of the NHIM and can now be referred to as a global TS trajectory in contrast to the local
TS trajectories as-sociated with the respective single barriers. The globalTS trajectory oscillates between the two local TS trajec-tories with the same frequency as the potential and canbe used to attach a mostly recrossing-free DS.The second set of central results of this work involvesthe determination of the decay rate constants of the os-cillatory (or time-dependent) two-saddle model system.Using the DS acquired in the first part, we can propagateensembles of particles, record a time-dependent reactantpopulation, and finally derive an instantaneous reactionrate parameterized by time according to Ref. [54]. Al-ternatively, the same result can be achieved purely byanalyzing the time-dependent phase space geometry. Forcomparison, a rate constant can be obtained from theglobal TS trajectory by means of Floquet stability anal-ysis. This method is in excellent agreement with theaverage of each instantaneous rate.While these results mark an important step in thetreatment of time-dependent multi-saddle systems, manyquestions still remain unanswered: First, we restrictedourselves to two-saddle systems with one DoF. To beapplicable to real-world systems, however, the methodspresented here will have to be generalized to at least moreDoFs because few chemical reactions can be treated ex-actly, or nearly so, when reduced to just one coordinate.Second, it will be important to investigate the influenceof minor manifold crossings on the rate constant. This isparticularly necessary for cases of time-dependent barri-ers found here in which there is no longer an equivalentto the primary geometric cross. This is even more chal-lenging when the alternation between barriers is drivenat high frequencies (cf. Fig. 2).The applicability of our results to real-world systemsalso remains to be demonstrated. It remains unclearwhether it is possible to treat systems without a pri-mary geometric cross. That chemical reactions can berepresented with potentials exhibiting the challenges dis-cussed here, is illustrated by the isomerization of ketenevia formylmethylene and oxirene which has been modeledvia a four-saddle potential [32, 58–61]. The isomerization0reaction of triangular KCN via a metastable linear K-CNconfiguration can similarly be described by a two-saddlesystem [62, 63]. Thus the analysis of time-dependentdriven potentials resolved here, when applied to theseand other chemical reactions, should provide new predic-tive rates for driven chemical reactions of interest.
ACKNOWLEDGMENTS
Useful discussions with Robin Bardakcioglu are grate-fully acknowledged. The German portion of thiscollaborative work was partially supported by theDeutsche Forschungsgemeinschaft (DFG) through GrantNo. MA1639/14-1. The US portion was partiallysupported by the National Science Foundation (NSF)through Grant No. CHE 1700749. MF is grateful for sup-port from the Landesgraduiertenf¨orderung of the LandBaden-W¨urttemberg. This collaboration has also ben-efited from support by the European Union’s Horizon2020 Research and Innovation Program under the MarieSk(cid:32)lodowska-Curie Grant Agreement No. 734557.
Appendix A: Tracking geometric crosses
To be able to track geometric crosses in phase spacereliably, the binary contraction method (BCM) (cf.Sec. II C) needs to be initialized with four points in fourdifferent reactive regions without human intervention.Since the geometric crosses of interest can be quite dis-torted while moving between saddles, the BCM cannotbe used precisely in the way described in Ref. [53]. Orig-inally, the initial quadrangle was defined by guessing thegeometric cross’s position and choosing two coordinateseach on a horizontal and a vertical line through it. Inour case, though, we first define an ellipse enclosing thegeometric cross, centered on this guess. The initial cor-ners of the quadrangle are then selected on the ellipseaccording to the algorithm described in Fig. 8.Interestingly, even a given geometric cross is not en-tirely free of substructures. Instead, it features a fractal-like set of crossing manifolds in its close proximity. Sincethese additional structures are extremely thin, however,they do not hinder the BCM from finding the desired geo-metric cross coordinates up to the desired precision whichmaybe be even smaller than the width of substructures.
Appendix B: Decay rates methods1. Ensemble method
The conceptually simplest method for calculating de-cay rates k e is by means of propagation of an ensemble.In analogy to Ref. [54] we first identify a line segmentparallel to the unstable manifold that satisfies the prop-erty: it lies on the reactant side between the stable man- π π π π π FIG. 8. Sketch of the algorithm for finding the inputs to theBCM: four arbitrary initial phase space coordinates repre-sentative of the four regions associated with the NHIM. In-stead of naively choosing coordinates on two right-angled axes(empty circles), an ellipse enclosing the geometric cross is con-structed. The first point (filled circle on right hand side) isselected arbitrarily at coordinates with angle α = 0. Sub-sequent coordinates are found by incrementing α in steps of∆ α = π/ α is temporarilyreduced until a point in the region in between is found. ifold and the DS at a distance that is small enough toallow for linear response and large enough to suppressnumerical instability. At t = t , an ensemble of particlesis placed on this line and propagated in time to yield atime-dependent reactant population N R ( t ; t ).In the second major step, one can obtain a reaction rateconstant k e by fitting an exponential decay N R ( t ; t ) ∝ exp[ − k e ( t − t )] to the reactant population [28, 64–68].This, however, is not possible in all systems—and, inparticular, in the two-saddle system being investigatedhere—because the decay in N R ( t ; t ) can be nonexponen-tial. Instead, we use the more general approach describedin Ref. [54], which involves examining the instantaneousdecays k e ( t ; t ) = − N R ( t ; t ) dd t N R ( t ; t ) . (B1)An analogous definition has been used in Ref. [14] forescape rates over a potential barrier.
2. Local manifold analysis
In general, the approach described above is compu-tationally expensive since a lot of particles have to bepropagated. This led us to develop a second method,called the LMA, for obtaining instantaneous reactionrates purely from the geometry of the stable and unstablemanifolds in phase space. If the slopes of the stable andunstable manifolds W s and W u at time t are given by∆ v s x ( t ) / ∆ x and ∆ v u x ( t ) / ∆ x then the instantaneous rate1 k m ( t ) can be written as their difference [54], k m ( t ) = ∆ v u x ( t ) − ∆ v s x ( t )∆ x . (B2)This allows for the calculation of k m ( t ) independently atdifferent times t , making it easy to compute in parallel.
3. Floquet method
Last, average decay rates k F can also be obtained di-rectly using a Floquet stability analysis [54, 56] that wasseen earlier to lead to accurate rates in reactions with atime-dependent barrier. While this method is the compu-tationally cheapest, it cannot yield instantaneous rates.To obtain the time-independent rate constant k F fora given trajectory on the NHIM, we can linearize theequations of motion using the Jacobian J ( t ) = (cid:18) − d V ( x, t ) (cid:14) d x (cid:19) . (B3)By integrating the differential equationd σ ( t )d t = J ( t ) σ ( t ) with σ (0) = , (B4)one obtains the system’s fundamental matrix σ ( t ). Whenconsidering trajectories with period T , M = σ ( T ) iscalled the monodromy matrix. Its eigenvalues m u and m s , termed Floquet multipliers, can be used to deter-mine the Floquet rate constant k F = 1 T (ln | m u | − ln | m s | ) . (B5) Appendix C: Possible errors for the discontinuousDS in Fig. 5(b)
There are multiple ways in which a discontinuous DScan show errors for a trajectory: (i)
A nonreactive trajectory is classified as reactive (la-beled “N as R” in Fig. 9). This can happen when a reac-tant enters the central region between the saddles whilethe DS is located near the left saddle. If this happensshortly before the DS jumps to the right, the particle can get reflected at the right barrier and leave the cen-tral region on the reactant side. (ii)
A reactive trajectory is classified as nonreactive(labeled “R as N” in Fig. 9). This can happen when areactant enters the central region while the DS is locatednear the right saddle. The DS can then jump discon-tinuously over the particle to the left saddle—which isnot counted as a crossing—before the particle leaves thecentral region to the right as a product. (iii)
The classification is thus inconsistent (labeled “in-cons.” in Fig. 9). This can happen, e. g., when a reactantenters the central region while the DS is located near the t p a r t i c l e f r a c t i o n [ % ] succ.N as RR as Nincons.recross. FIG. 9. Fraction of trajectories with (see text) and without(labeled “succ.”) errors as a function of time t for the dis-continuous DS from Fig 5. right saddle, and leaves it again to the reactant side whilethe DS is near the left saddle. As a result, a backwardreaction is recorded even though the particle started asa reactant. Similarly, it is possible to detect two for-ward reactions over the same DS without an intermedi-ate backward reaction. It is unclear which reaction is tobe counted as the real one. (iv) The trajectory crosses the DS multiple times, i. e.,it exhibits recrossings (labeled “recross.” in Fig. 9). Thiscan happen when a particle enters and leaves the centralregion on the reactant side while the DS is near the leftsaddle, resulting in two crossings.The distributions of the errors that arise from the dis-continuous DS originating from recrossing and misclas-sification are shown in Fig. 9. The sum of these errorsgives rise to that reported in Fig. 5 for this DS. [1] D. G. Truhlar, W. L. Hase, and J. T. Hynes, J. Phys.Chem. , 2664 (1983), doi:10.1021/j100238a003.[2] R. A. Marcus, Science , 1523 (1992),doi:10.1126/science.256.5063.1523.[3] D. G. Truhlar, B. C. Garrett, and S. J. Klippenstein, J. Phys. Chem. , 12771 (1996), doi:10.1021/jp953748q.[4] T. Komatsuzaki and R. S. Berry, Proc. Natl. Acad. Sci.U.S.A. , 7666 (2001), doi:10.1073/pnas.131627698.[5] E. V. Anslyn and D. A. Dougherty, Modern Physical Or-ganic Chemistry (University Science Books, Mill Valley, CA, 2006).[6] M. R. Harper, K. M. Van Geem, S. P. Pyl, G. B. Marin,and W. H. Green, Combustion and Flame , 16 (2011),doi:10.1016/j.combustflame.2010.06.002.[7] S. T. Chill, J. Stevenson, V. Ruehle, C. Shang, P. Xiao,J. D. Farrell, D. J. Wales, and G. Henkelman, Journalof Chemical Theory and Computation , 5476 (2014),doi:10.1021/ct5008718.[8] Z. C. Kramer, B. K. Carpenter, G. S. Ezra, andS. Wiggins, J. Phys. Chem. A , 6611 (2015),doi:10.1021/acs.jpca.5b02834.[9] B. Chen, R. Hoffmann, and R. Cammi, Angew. Chem.,Int. Ed. , 11126 (2017), doi:10.1002/anie.201705427.[10] W. H. Miller, J. Chem. Phys. , 1651 (1968),doi:10.1063/1.1668891.[11] W. H. Miller, J. Phys. Chem. , 960 (1979),doi:10.1021/j100471a015.[12] C. R. Doering and J. C. Gadoua, Phys. Rev. Lett. ,2318 (1992), doi:10.1103/PhysRevLett.69.2318.[13] J. Maddox, Nature (London) , 771 (1992),doi:10.1038/359771a0.[14] J. Lehmann, P. Reimann, and P. H¨anggi, Phys. Rev.Lett. , 1639 (2000), doi:10.1103/PhysRevLett.84.1639.[15] J. Lehmann, P. Reimann, and P. H¨anggi, Phys. Rev. E , 6282 (2000), doi:10.1103/PhysRevE.62.6282.[16] J. Lehmann, P. Reimann, and P. H¨anggi, Phys. StatusSolidi B , 53 (2003), doi:10.1002/pssb.200301774.[17] C. Van den Broeck, Phys. Rev. E , 4579 (1993),doi:10.1103/PhysRevE.47.4579.[18] P. H¨anggi, Chem. Phys. , 157 (1994),doi:10.1016/0301-0104(93)E0422-R.[19] T. D. Shepherd and R. Hernandez, J. Chem. Phys. ,2430 (2001), doi:10.1063/1.1386422.[20] T. D. Shepherd and R. Hernandez, J. Chem. Phys. ,9227 (2002), doi:10.1063/1.1516590.[21] T. D. Shepherd and R. Hernandez, J. Phys. Chem. B , 8176 (2002), doi:10.1021/jp020620h.[22] J. M. Moix, T. D. Shepherd, and R. Hernandez, J. Phys.Chem. B , 19476 (2004), doi:10.1021/jp046629w.[23] R. D. Vogelaere and M. Boudart, J. Chem. Phys. ,1236 (1955), doi:10.1063/1.1742248.[24] W. H. Miller, J. Chem. Phys. , 2216 (1976),doi:10.1063/1.433379.[25] E. Pollak, M. S. Child, and P. Pechukas, J. Chem. Phys. , 1669 (1980), doi:10.1063/1.439276.[26] B. K. Carpenter, “Potential energy surfaces and reactiondynamics,” in Reactive Intermediate Chemistry (JohnWiley & Sons, New York, 2005) Chap. 21, pp. 925–960.[27] P. Pechukas,
Dynamics of Molecular Collisions , edited byW. H. Miller, Vol. 2 (Plenum, New York, 1976) Chap. 6,pp. 269–322.[28] P. Pechukas, Annu. Rev. Phys. Chem. , 159 (1981),doi:10.1146/annurev.pc.32.100181.001111.[29] E. Pollak and P. Pechukas, J. Chem. Phys. , 325(1979), doi:10.1063/1.437194.[30] J. Rehbein and B. K. Carpenter, Phys. Chem. Chem.Phys. , 20906 (2011), doi:10.1039/C1CP22565K.[31] P. Collins, B. K. Carpenter, G. S. Ezra, andS. Wiggins, J. Chem. Phys. , 154108 (2013),doi:10.1063/1.4825155.[32] G. T. Craven and R. Hernandez, Phys. Chem. Chem.Phys. , 4008 (2016), doi:10.1039/c5cp06624g.[33] A. Junginger, L. Duvenbeck, M. Feldmaier, J. Main,G. Wunner, and R. Hernandez, J. Chem. Phys. , 064101 (2017), doi:10.1063/1.4997379.[34] S. Glasstone, K. J. Laidler, and H. Eyring, The The-ory of Rate Processes: The Knetics of Chemical Reac-tions, Viscosity, Diffusion and Electrochemical Phenom-ena (McGraw-Hill, New York, 1941).[35] K. Fukui, J. Phys. Chem. , 4161 (1970),doi:10.1021/j100717a029.[36] S. Kato and K. Fukui, J. Am. Chem. Soc. , 6395(1976), doi:https://doi.org/10.1021/ja00436a061.[37] H. Eyring, J. Chem. Phys. , 107 (1935),doi:10.1063/1.1749604.[38] E. P. Wigner, J. Chem. Phys. , 720 (1937),doi:10.1063/1.1750107.[39] M. Feldmaier, P. Schraft, R. Bardakcioglu, J. Reiff,M. Lober, M. Tsch¨ope, A. Junginger, J. Main,T. Bartsch, and R. Hernandez, J. Phys. Chem. B ,2070 (2019), doi:10.1021/acs.jpcb.8b10541.[40] A. J. Lichtenberg and M. A. Liebermann, Regular andStochastic Motion (Springer, New York, 1982).[41] R. Hernandez and W. H. Miller, Chem. Phys. Lett. ,129 (1993), doi:10.1016/0009-2614(93)90071-8.[42] R. Hernandez, Ph.D. thesis, University of California,Berkeley, CA (1993).[43] E. Ott,
Chaos in Dynamical Systems , 2nd ed. (CambridgeUniversity Press, Cambridge, England, 2002).[44] S. Wiggins,
Normally Hyperbolic Invariant Manifolds inDynamical Systems , Vol. 105 (Springer Science & Busi-ness Media, New York, 2013).[45] T. Komatsuzaki and R. S. Berry, Phys. Chem. Chem.Phys. , 1387 (1999), doi:10.1039/A809424A.[46] C.-B. Li, A. Shoujiguchi, M. Toda, and T. Komatsuzaki,Few-Body Systems , 173 (2006).[47] C.-B. Li, A. Shoujiguchi, M. Toda, and T. Ko-matsuzaki, Phys. Rev. Lett. , 028302 (2006),doi:10.1103/PhysRevLett.97.028302.[48] C.-B. Li, M. Toda, and T. Komatsuzaki, J. Chem. Phys. , 124116 (2009), doi:10.1063/1.3079819.[49] C. Mendoza and A. M. Mancho, Phys. Rev. Lett. ,038501 (2010), doi:10.1103/PhysRevLett.105.038501.[50] A. M. Mancho, S. Wiggins, J. Curbelo, and C. Mendoza,Commun. Nonlinear Sci. Numer. Simul. , 3530 (2013),doi:10.1016/j.cnsns.2013.05.002.[51] G. T. Craven and R. Hernandez, Phys. Rev. Lett. ,148301 (2015), doi:10.1103/PhysRevLett.115.148301.[52] A. Junginger and R. Hernandez, J. Phys. Chem. B ,1720 (2016), doi:10.1021/acs.jpcb.5b09003.[53] R. Bardakcioglu, A. Junginger, M. Feldmaier, J. Main,and R. Hernandez, Phys. Rev. E , 032204 (2018),doi:10.1103/PhysRevE.98.032204.[54] M. Feldmaier, R. Bardakcioglu, J. Reiff, J. Main, andR. Hernandez, J. Chem. Phys. , 244108 (2019),doi:10.1063/1.5127539.[55] M. Feldmaier, J. Reiff, R. M. Benito, F. Borondo,J. Main, and R. Hernandez, J. Chem. Phys. , 084115(2020), doi:10.1063/5.0015509.[56] G. T. Craven, T. Bartsch, and R. Hernandez, J. Chem.Phys. , 041106 (2014), doi:10.1063/1.4891471.[57] G. S. Hammond, J. Am. Chem. Soc. , 334 (1955),doi:10.1021/ja01607a027.[58] J. D. Gezelter and W. H. Miller, J. Chem. Phys. ,7868 (1995), doi:10.1063/1.470204.[59] I. S. Ulusoy, J. F. Stanton, and R. Hernandez, J. Phys.Chem. A , 7553 (2013), doi:10.1021/jp402322h.[60] F. A. L. Maugi`ere, P. Collins, G. Ezra, S. C. Farantos, and S. Wiggins, Theor. Chem. Acta , 1507 (2014),doi:10.1007/s00214-014-1507-4.[61] I. S. Ulusoy and R. Hernandez, Theor. Chem. Acc. ,1528 (2014), doi:10.1007/s00214-014-1528-z.[62] H. P´arraga, F. J. Arranz, R. M. Benito, andF. Borondo, J. Chem. Phys. , 194304 (2013),doi:10.1063/1.4830102.[63] H. P´arraga, F. J. Arranz, R. M. Benito, andF. Borondo, J. Phys. Chem. A , 3433 (2018),doi:10.1021/acs.jpca.8b00113.[64] D. Chandler, J. Chem. Phys. , 2959 (1978),doi:10.1063/1.436049. [65] W. H. Miller, S. D. Schwartz, and J. W. Tromp, J. Chem.Phys. , 4889 (1983), doi:10.1063/1.445581.[66] P. H¨anggi, P. Talkner, and M. Borkovec, Rev. Mod.Phys. , 251 (1990), doi:10.1103/RevModPhys.62.251,and references therein.[67] K. A. Connors, Chemical Kinetics: The Study of Reac-tion Rates in Solution (Wiley-VCH, Weinham, Germany,1990).[68] S. K. Upadhyay,