Thermal decay rates of an activated complex in a driven model chemical reaction
Robin Bardakcioglu, Johannes Reiff, Matthias Feldmaier, Jörg Main, Rigoberto Hernandez
TThermal decay rates of an activated complex in a driven model chemical reaction
Robin Bardakcioglu, 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: January 6, 2021)Recent work has shown that in a non-thermal, multidimensional system, the trajectories in theactivated complex possess different instantaneous and time-averaged reactant decay rates. Underdissipative dynamics, it is known that these trajectories, which are bound on the normally hyper-bolic invariant manifold (NHIM), converge to a single trajectory over time. By subjecting thesedissipative systems to thermal noise, we find fluctuations in the saddle-bound trajectories and theirinstantaneous decay rates. Averaging over these instantaneous rates results in the decay rate of theactivated complex in a thermal system. We find, that the temperature dependence of the activatedcomplex decay in a thermal system can be linked to the distribution of the phase space resolveddecay rates on the NHIM in the non-dissipative case. By adjusting the external driving of the reac-tion, we show that it is possible to influence how the decay rate of the activated complex changeswith rising temperature.
I. INTRODUCTION
In order to predict the rate of chemical reactions,transition state theory (TST) [1–15] utilizes a dividingsurface (DS) in phase space [16, 17] to determine when re-actants decay into products under quasi-equilibrium con-ditions [18]. We and others [10–15, 19–24] have extendedthe use of TST to address reactions far out of equilibriumleading to rates, resolution of the reactive geometry, orthe reaction paths. Such conditions can arise when thereactants respond to external stimuli—e. g. under drivenconditions or collective effects of the reacting environ-ment.The accuracy of the TST rate depends on the accuracyof the DS to truly divide the space between reactantsand products. That is, it must satisfy the condition thatno reacting particle crosses it more than once. Unfor-tunately, in a driven system, it is often not enough toclassify molecules based only by their spatial configura-tion in large part because the structure of the reactiongeometry is itself time-dependent. The DS must thenbe extended to the full phase space of the system in atime-dependent frame to ensure that it is free of recross-ings [19, 20, 25, 26].An additional complexity arises in activated processeswhich must either explicitly address a solvent by extend-ing the number of degrees to a macroscopic degree—e. g.moles—or implicitly by including them through a for-malism such as that of Brownian motion [27–29]. Ineither case, the coupling between the reactants and thesolvent can be surmised through an effective friction. Ac-tivated processes have been seen to undergo a Kramersturnover [28, 29] in the rate as they are solvated fromlow to high friction [30–34], and hence the exact value ∗ Correspondence to: [email protected] of the friction is important in determining the rate. Ingeneral, the activated complex is a collection of unsta-ble configurations located near the energy barrier be-tween the reactants and products. In the modern lan-guage of differential geometry, this has become associatedwith the so-called normally hyperbolic invariant mani-fold (NHIM) [35–40]. It is a co-dimension 2 manifold inthe phase space of the system, which is characterized bythe condition that trajectories started on that manifoldstay there when propagated forward or backward in time.The DS (or NHIM) is complementary to the brute forcecalculation of reaction rates using trajectories that startfrom the reactant region. Such trajectories only countif they cross to products and are rare in activated pro-cesses. Nevertheless, the few that are reactive, cross anexact DS once and only once. Avoiding the work of deter-mining the nonreactive trajectories from the reactant re-gion, one thus typically focuses on the rate of trajectoriesleaving from the DS which when exact—because of thenon-recrossing condition—gives the flux-over-populationrate [41–44]. Either because of time dependence or di-mensionality, the NHIM itself can generally accommo-date a set of trajectories that neither enter nor leave it.In recent work, we have explored the stability of thisclass of trajectories as we have conjectured that theirdecay is connected to the decay of the reactive trajec-tories [20, 40, 45]. In dissipative driven systems, due tothe properties of the NHIM, trajectories within it con-verge towards a single transition state (TS) trajectoryafter a sufficiently long time in the saddle region. Thosetrajectories near it, will also be trapped towards a singletrajectory of the NHIM [12, 23, 46]. Obtaining the decayrate of the trajectories within the NHIM, and a determi-nation of how the thermal environment affects them isthe primary contribution of this work.Using the system and methods described in Sec. II,we can use geometric structure obtained directly todetermine rates in driven chemical reactions that are a r X i v : . [ phy s i c s . c h e m - ph ] J a n not isolated, but rather coupled to a dissipative en-vironment. This is a necessary advance for the useof the nonrecrossing dividing surfaces— viz. the time-dependent DS attached to the NHIM—that we and oth-ers [12, 20, 23, 40, 45–47] have been developing becausemany chemical reactions of interest occur in a solvent.The results presented in Sec. III provide a demonstrationof the stochastic time-dependent motion of the NHIM atfixed orthogonal modes (Sec. III A) and the collapse ofthe transition states under dissipation towards a singletrajectory on the NHIM (Sec. III B). The time depen-dence of the instantaneous reactant decay rate and thetemperature dependence of the average decay rate of theactivated complex over long times are shown in Sec. III C.The temperature-dependent behavior of the average de-cay rate is linked to the phase space resolved average de-cay rate of the non-thermal system. We also find that thetemperature dependence of the decay rate can be influ-enced by changing the oscillation of the periodic driving. II. SYSTEM AND METHODS
In this section, we first recapitulate the representationof a chemical reaction [20, 40, 48]. As in earlier work,we impose Langevin dynamics to represent the influenceof the bath, and use a driven saddle potential to revealthe decay rates of trajectories within the NHIM. Thespecifics of the system and the associated equation ofmotion (EOM) are summarized in Sec. II A. The unsta-ble transition states, i. e., the trajectories on the NHIM,can then be constructed using the approach described inSec. II C [45, 49]. The instantaneous decay rates at coor-dinates of the NHIM and the average decay rates along atransition state can then be constructed as summarizedin Sec. II D.
A. Model chemical reaction
The dynamics of the system under investigation isgiven by the Langevin equation,˙ v = − γ v + ξ ( t ) − ∇ V ( x , t ) , (1a)˙ x = v , (1b)where x is the coordinate vector of the system, v thecorresponding velocities, t the time, and γ the frictioncoefficient. The vector ξ ( t ) represents the fluctuationsaround the time-dependent mean-force potential V ( x , t ).Here, we represent each component as white noise whichsatisfies the fluctuation-dissipation theorem [50–52] withrespect to the specified friction, h ξ i ( t ) i = 0 , (2) h ξ i ( t ) ξ j ( t ) i = 2 γT δ ij δ ( t − t ) , (3)where δ ij and δ ( t − t ) represent the Kronecker delta andDirac delta distribution, respectively. The temperature T is given in units where the Boltzmann constant k B is1 so as to give it the same units as energy. Note that thesame noise sequence is used for all trajectories.The specific potential investigated here has been usedmany times before in previous work [19, 20, 40, 45, 48, 49,53, 54]. It is a two-dimensional rank-1 saddle potentialof the form V ( x, y, t ) = 2 exp (cid:16) − [ x − . ωt )] (cid:17) + 2 (cid:18) y − π arctan (2 x ) (cid:19) . (4)This potential models a reaction over an energy barrieroscillating with frequency ω and provides a nonlinear cou-pling between reaction coordinate and orthogonal modealong the reaction path. The coordinate component x ap-proximates the reaction coordinate, i. e., the unstable di-rection of the saddle potential, by construction. Likewise,the y component approximates the orthogonal modes.These coordinates x = ( x, y ), together with their veloci-ties ˙ x = v x , ˙ y = v y form the phase space of the system.The effect of the Langevin dynamics on the trajectoriesof the system is illustrated in Fig. 1. The phase spacecoordinates of two arbitrary TS trajectories on the NHIMare shown. In the non-thermal case we can see that thetrajectory follows a smooth path in phase space and itis periodic in sync with the deterministic driving. Incontrast, the trajectory of the thermal system fluctuatessignificantly, especially in the velocities. The stochasticand aperiodic motion of the thermal trajectories presentsa new challenge for our earlier methods on non-thermalsystems [20, 40, 45, 49] addressed here. B. Identification of the NHIM
The NHIM, barring any general definition and hereonly limited to a rank-1 saddle potential, is the set of alltrajectories in phase space bound forever to the saddleregion in both forward and backward in time directions.It is a manifold of co-dimension two in phase space. It isalso associated with a pair of stable and unstable mani-folds of co-dimension one, whose closures intersect at theNHIM. For our model chemical reaction [Eq. (4)], theposition ( x, v x ) NHIM ( y, v y ) of the NHIM can be param-eterized as a function of the stable orthogonal modes y and v y at a specific time t .The motion of individual trajectories in a close neigh-borhood of the NHIM is stochastic for a thermal systemas illustrated in Fig. 1. Nevertheless, the correspond-ing stable and unstable manifolds in Fig. 2(a) remainsmooth, and generally so. Thus, even for a thermalsystem, the phase space in a local neighborhood lookssimilar to a non-thermal system, as also observed inRefs. [20, 40, 45, 49]. Hence, the typical cross-like inter-section of the stable and unstable manifolds is preservedand the position of the NHIM can be directly obtained − . . x γ = 0 . T = 0 . γ = 0 . T = 1 . − . . y − v x . . . . . t − v y FIG. 1. Phase space coordinates ( x, y, v x , v y ) of two trajec-tories arbitrarily selected from the ensemble of TS trajecto-ries on the NHIM at T = 0 and T = 1 .
0, respectively. Thenon-thermal trajectory ( T = 0, γ = 0) is periodic with thedriving frequency ω = π of the potential, as can be readilyseen. The thermal trajectory ( T = 1 . γ = 0 . from the intersection of the stable and the unstable man-ifolds. In a thermal system, however, the intersection( x, v x ) NHIM ( y, v y ) will not only depend on the specificchoice of the orthogonal modes ( y, v y ), but also on the pa-rameters γ and T . This finding is illustrated in Fig. 2(a)for the paradigmatic system used throughout this work.The stable and unstable manifolds have the char-acteristic property that any trajectory near them ispropagated towards or away from the NHIM, respec-tively. Only at the intersection of their closures— viz. theNHIM— trajectories are unstably bound to the saddleregion. As the stable and unstable manifold are them-selves invariant manifolds, they cannot be crossed by anytrajectory of the system. This effectively separates thephase space in the close neighborhood of the NHIM into − . − . − . . . . . x − − v x (a) γ = 0 . , T = 0 .
0, stable γ = 0 . , T = 0 .
0, unstable γ = 0 . , T = 1 .
0, stable γ = 0 . , T = 1 .
0, unstable (b) v x x I IV IIIII
FIG. 2. (a) Stable and unstable manifold at orthogonal modes y = v y = 0 for a non-thermal and thermal system with driv-ing frequency ω = π . The intersection of these manifoldsmark the x and v x coordinate of the NHIM for these orthog-onal modes. Crosses such as these can be obtained for anyorthogonal modes ( y, v y ) and time t . Thus, one can expressthe reaction coordinates ( x, v x ) NHIM ( y, v y , t ) of the NHIM asfunctions of initial time and orthogonal modes. (b) Sketch ofthe stable and unstable manifold. The arrowheads indicatethe projected path trajectories take within the four separatedregions. A quadrangle with a vertex in each of the four re-gions illustrates the binary contraction method, an iterativeroutine to find the NHIM at the intersection of the manifolds,see Sec. II B. four distinct regions demarcated by these manifolds, seeFig. 2(b). For a chosen saddle region x ∈ [ x min , x max ]with appropriately chosen boundaries x min and x max , thedynamics of trajectories is characteristic for any of theseregions, as indicated with the black arrows in Fig. 2(b).In region (I), non-reactive trajectories originate from thereactant side x < x min in the past and fall back to thereactant side in the future. The same is true for region(III) which contains all trajectories that originate fromthe product side x > x max and also fall back to the prod-uct side. Regions (II) and (IV) hold the reactive trajecto-ries from the reactant to the product side or, respectively,vice versa.We use the binary contraction method (BCM) [49] tonumerically find the NHIM at a given time t for a givenset ( y, v y ) of orthogonal modes. The procedure of theBCM is initiated with a quadrangle having its four ver-tices (blue dots in Fig. 2(b)) in each of the four reac-tive and non-reactive regions in a close neighborhood ofthe NHIM. In successive steps, the midpoint betweenpairs of nearby points of the quadrangle is propagatedto determine which region it belongs to and then re-places the corresponding point. The area of the quad-rangle thus shrinks and its center converges to the inter-section ( x, v x ) NHIM ( y, v y ) of the manifold. The resulting( y, v y ) corresponds to the position ( x, v x ) of the NHIM.This convergence is exponentially fast and, therefore, theBCM is very efficient in finding the NHIM for a given setof orthogonal modes. For the technical details of theBCM, we refer the reader to Ref. [49]. C. Trajectories on the NHIM
Due to the hyperbolic nature of the NHIM, trajecto-ries on the NHIM are unstable. Small deviations fromthe NHIM will grow exponentially in time until the tra-jectory leaves its immediate vicinity. Thus any deviationin a point relative to the NHIM, no matter how small,will lead its subsequent propagation to fall off of it even-tually. This presents a challenge to the propagation oftrajectories on the NHIM using numerical simulations.These numerical deviations can be suppressed througha machine-learning representation of the NHIM project-ing them back onto the NHIM as has been done inRef. [45]. Here, we employ this approach leveraging thenumerically accurate BCM presented in Ref. [49]. Giventhe system and bath parameters, the reaction coordinateand corresponding velocity of the NHIM can be inte-grated in time as a function of the remaining orthogo-nal modes in the usual way. By subsequently project-ing the trajectory back onto the NHIM it is effectivelypropagated in a system with reduced dimensions whichis spanned by the orthogonal modes. The result is a tra-jectory that remains on the NHIM.
D. Decay rates of reactant population
The relative stability of the NHIM can be character-ized through the instantaneous decay rates k ( y, v y , t ) ofthe reactant population for the specific orthogonal modes( y, v y ) in its close neighborhood at each time t . Thesedecay rates can be obtained directly by propagating en-sembles of reactive trajectories but this can become cum-bersome and numerically expensive. The local manifoldanalysis (LMA) was introduced in Ref. [40] to overcomethis problem by leveraging the linear dynamics near theNHIM.Since the dynamics relative to the transition state islinear, it is possible to propagate trajectories using a lin-ear map, the fundamental matrix M , obtained via˙ M = J ( t ) M , (5) with the initial condition of M ( t ) = , and a Jacobian J ( t ) = ∂ ( ˙ x , ˙ v ) ∂ ( x , v ) = − ∂ V∂x − ∂ V∂x∂y − γ − ∂ V∂x∂y − ∂ V∂y − γ (6)parameterized in time along a specific trajectory of thetransition state. The stochastic force ξ ( t ) does not con-tribute to any component of the Jacobian in Eq. (6) as itis purely time-dependent. However, this does not meanthat it is neglected in the Langevin dynamics. The Ja-cobian (6) describes the motion relative to a transitionstate whose trajectory is fluctuating under the influenceof the stochastic force. That means that although thedynamics relative to the close vicinity of the transitionstate might not be fluctuating, the total dynamics stilldoes. Using the linear dynamics near the transition stateaccording to Eq. (5), we can extract the instantaneousmotion of a particle ensemble near that state from theJacobian.The LMA models the presence of such a uniform par-ticle ensemble in the close neighborhood of the NHIM.For the system of Eqs. (1) and (4) at specific orthogonalmodes ( y, v y ) and time t , we can obtain a local instanta-neous decay rate k ( y, v y , t ) = ∆ v u x ∆ x u ( y, v y , t ) − ∆ v s x ∆ x s ( y, v y , t ) , (7)where ∆ v s , u x / ∆ x s , u represents the slopes of the stable andunstable manifolds in the corresponding ( x, v x ) cross sec-tion of the full phase space near the transition state. De-termining these slopes, or rather the instantaneous decayrate as the difference of these slopes, is the main objec-tive of the numerical implementation of the LMA. Aderivation of the LMA can be found in Ref. [40], and theadditional corrections that would be necessary for moregeneral cases can be found in the Supplementary Mate-rial of Ref. [55].An average rate for the decay relative to trajectories onthe NHIM is obtained by computing the time average ofthe instantaneous decay rates for a sufficiently long timeto obtain convergence. For a trajectory that is initializedat time t and position ( y , v y ) on the NHIM, and pa-rameterized by the orthogonal modes y ( t ) and v y ( t ), thisaverage yields¯ k ( y , v y , t ) = lim τ →∞ τ Z t + τt k ( y ( t ) , v y ( t ) , t ) d t . (8)In the special case of T = 0, we find that the trajectoriesare periodic or quasi-periodic, and it suffices to integratefor the period or quasi-period. An alternative approachto obtain mean decay rates is provided by a Floquet anal-ysis of said trajectory [40, 56]¯ k ( y , v y , t ) = lim τ →∞ τ (ln | m l ( τ ) | − ln | m s ( τ ) | ) , (9) . . . . . t − . − . . x NH I M ( y , v y , t ) γ = 0 . T = 0 . γ = 0 . T = 1 . FIG. 3. Position x NHIM of the NHIM for fixed y = v y = 0and ω = π as a function of time. Without the Langevinterms, the NHIM (dashed curve) oscillates periodically, as onewould expect from the periodically driven saddle potential.With the Langevin terms, however, the NHIM (solid curve)moves stochastically. Compared to a trajectory on the NHIMfor the same parameters, the reaction coordinates x does notfluctuate as much as the NHIM in phase space, see Fig. 1. where m l , s ( t ) are the eigenvalues of the fundamental ma-trix M . It reduces to the monodromy matrix when thetrajectory is periodic. Here, the subscripts l and s denotethe eigenvalues with the largest and smallest absolute val-ues, respectively [40, 56]. III. RESULTS AND DISCUSSIONA. Stochastic motion of the NHIM under noise
The time-dependent NHIM and the associated decayrates in a driven chemical reaction are resolved here forthe model system of Eq. (1). If the dissipative andstochastic terms are excluded, the result is a smoothoscillating motion with the same period as that of theoscillating potential [55]. However, when the system issubject to the Langevin terms— viz. friction and thermalnoise—the motion of the NHIM becomes stochastic asillustrated in Fig. 3. Here, the expected behavior in thetime-dependence of x NHIM for a specific set of coordinates( y, v y ) in the dynamics without and with the Langevinterms is apparent. The latter case now exhibits stochasticfluctuations in both position (as shown) and momentum(as not shown) space. They are in response to the com-bination of the collective stochastic thermal driving andthe periodic driving terms. Such fluctuations were notseen in the position space of the TS trajectories shown inFig. 1 because in this relatively weak friction regime, thehigh-frequency fluctuations are very small. However, inFig. 3, the NHIM does exhibit short-time fluctuations inthe position space as a manifestation of the overall phasespace motion. Nevertheless, the DS is recrossing-free asit incorporates the noise. − − x − − y (a) − − x − − y (b) FIG. 4. The PSOSs of the system defined by Eqs. (1) and(4) for varying friction at temperature T = 0 and saddle fre-quency ω = π . The non-thermal system with friction γ = 0and the thermal system with γ = 0 . B. Dissipative dynamics on the NHIM
Although we have seen that the dynamics of the systemin the general case of Eq. (1) for the potential in (4) isnot periodic, it is nevertheless instructive to examine thestroboscopic Poincar´e surface of section (PSOS) of its tra-jectories. Similar to the approach used in Refs. [40, 48],we record the position of trajectories in the ( y, v y ) sectionin phase space at time steps equal to the period of thedriving. As we start with a point on the NHIM, it neces-sarily must stay on the NHIM. Thus, the coordinates inthe PSOS, as seen e. g. in Fig. 4, remain projected ontothe NHIM.The PSOSs of Fig. 4 show the contrast in the dynamicsupon the introduction of friction. In the upper panel,the dynamics on the NHIM is regular. It gives rise tothe expected stable concentric tori and an elliptic fixedpoint. In the bottom panel, as a result of the friction, t k ( t ) (a) T = 0 . T = 0 . t k ( t ) (b) FIG. 5. (a) Comparison of instantaneous decay rates of theequilibrium trajectory on the NHIM as defined in Sec. III C forfriction γ = 0 . T = 0 . ω = π . The instanta-neous rate of the thermal complex roughly follows the oscilla-tions of the non-thermal trajectory but exhibits fluctuationsto a certain degree. (b) Instantaneous rate of the thermalactivated complex in over a longer time interval. Stochasticfluctuations of the rate become more evident over hundredsof saddle oscillation periods. The dashed line highlights thetime-average of the rate. the would-be tori now spiral towards a fixed point inthe stroboscopic projection. This fixed point refers to atime-periodic trajectory in phase space, which acts as anattractor for particles on the NHIM. C. Thermal decay rates
The collapse of the NHIM to a single periodic trajec-tory in dissipative regimes with temperature T = 0 ob-served in Sec. III B can be used to our advantage fortemperatures above zero. Even when noise is in play, thefluctuating trajectories on the NHIM will approach a sin-gle fluctuating trajectory over long times, which we willrefer to as the equilibrium trajectory on the NHIM. It canbe determined numerically by propagating an arbitrarypoint on the NHIM for a sufficiently long time until theequilibrium is reached. The initial build-up is discardedwhen calculating rates. In the infinite time limit, use ofthe equilibrium trajectory to obtain the average decayrate ¯ k eliminates its dependence on the initial conditions( y , v y , t ) in Eqs. (8) and (9). That is, all contributions to the average decay rate that would depend on the ini-tial conditions are dwarfed by the contributions of theequilibrium trajectory.Assuming that the long-term behavior is independentof the initial time t , the equilibrium trajectory can beused to construct the expected value h k i ( T, γ ) of the re-actant decay rate as a function of the temperature T andfriction γ h k i ( T, γ ) = ¯ k ( y , v y , t ; T, γ ) , (10)which, due to the fact that the initial conditions will beforgotten by the dynamics over time, is the same for anyset of initial conditions ( y , v y , t ).
1. Instantaneous decay rates
The influence of the noise on the average decay rate isrevealed by the time evolution of the instantaneous de-cay rate. The rates shown in Fig. 5(a) are plotted overfive saddle oscillations. Despite the stochastic nature ofthe NHIM at fixed orthogonal modes ( y, v y ), we find thatthe instantaneous rate of the thermal trajectory on theNHIM is smooth and still roughly follows the regular os-cillation we find for zero temperature. This can be at-tributed to the fact that the trajectories, as integratedvalues over a noisy acceleration, have a smooth time evo-lution. Despite the instant decay in the time correlationof the stochastic force in the fluctuation-dissipation rela-tion [Eq. (3)], these trajectories have a finite memory oftheir previous positions [57]. However, over sufficientlylong time scales, it is possible to obtain a time evolutionof the instantaneous reactant decay rate that resemblesan uncorrelated fluctuation, as can be seen in Fig. 5(b).
2. Temperature dependence of average decay rates
The averaged rate has to be determined over severalhundreds—sometimes even a few thousand—saddle oscil-lations to obtain a statistically-sound converged average.This is as a consequence of the need to achieve the lim-iting equilibrium trajectory as observed in Sec. III C 1.The resulting temperature dependence of the averageddecay rate for two sets of parameters is shown in Fig. 6.These parameters were chosen, as they give rise to twovery different regimes in the shape of the rate from con-cave to convex. With a driving frequency of ω = 0 . π ,we obtain a rate that is monotonically rising, as temper-ature rises. However, at a driving frequency of ω = π , weobtain a rate that at first decreases, before it increaseswith rising temperature.The change in behavior of the rate curves may seemcounter-intuitive at first. One might expect that a highertemperature—i. e., a higher average energy—would causethe reaction to surmount an energy barrier at a higherrate. This expectation is not contradicted by the present . . . . T . . . h k i ω = 1 . πω = 0 . π FIG. 6. Temperature dependent time-average of the instan-taneous decay rate for the activated complex of the modelsystem at two different driving frequencies, ω = 0 . π and1 . π . The dotted vertical line highlights the average rates at T = 0 . results. The decay rate computed here is the decay rateof the reactant population close to the transition state,i. e., the rate of reaction after the reactant has alreadysurmounted the energy barrier. Such a rate does not in-clude the increased population of activated reactants thatwould arise from a higher temperature and consequentlydoes not have to increase accordingly.We find that strong noise causes the activated complexto fluctuate strongly near the NHIM as seen in Fig. 7.When comparing these coordinate fluctuations with thephase space resolved decay rates of the non-thermal sys-tem according to Fig. 7 or additional examples in theSupplemental Material [58], we can also see how a ther-mal activated complex explores several transition statesof the non-thermal system. This further suggests thatthe decay of the thermal activated complex into reactantsor products is related to the distribution of the explored,non-thermal transition states. This interpretation is con-sistent with the observation made in Sec. II D that noiseonly affects the average decay rate by the trajectory thatis used to obtain it. Even though this analogy is notformally exact since the Jacobian J according to Eq. (6)contains the friction coefficient γ —which is not zero inthe thermal case—the conjectured interpretation appearsto hold for the low-temperature thermal case discussedin Figs. 6 and 7.A heuristic argument in support of the conjecture is asfollows. As temperature rises, the activated complex de-viates further from the corresponding periodic trajectoryat T = 0. This in turn causes the activated complex toexplore transition states that are further from said tra-jectory. Moreover, the distribution of explored transitionstates expands into a region with shrinking decay ratesas temperature rises. The equilibrium trajectory now re-sides in a local maximum of decay rates, as is the case for ω = 0 . π , and hence the decay of the activated complexdecreases with rising temperature. − . − . . . . y − . − . − . . . v y (a) . . . h k i − . − . . . . y − . . . . . v y (b) . . . . h k i FIG. 7. Stroboscopic projections of a typical trajectory on theNHIM for each of the saddle driving frequencies, (a) ω = π and (b) ω = 0 . π , in the dissipated thermal case, γ = 0 . T = 0 .
2, are shown as white filled circles. Refer to Fig. 6 forthe corresponding average rates where the vertical dashed linecrosses the curves. The color maps show the time-averageddecay rates h k i of a non-thermal activated complex on theNHIM for the initial conditions ( y, v y ) at time t = 0. Forreference, the fixed points of the stroboscopic maps in thenon-thermal case are indicated by a cross. IV. SUMMARY AND CONCLUSION
We have characterized the geometric structure of amodel chemical reaction, thereby taking into accountboth external driving and noise and friction describedby the Langevin terms. We have shown that the tem-perature dependence of the activated complex decay inthis thermal system is linked to the distribution of thephase space resolved decay rates on the NHIM in thenon-dissipative case. The decay rate of the activatedcomplex depends on the external driving and the tem-perature, and these dependencies can be used to controlthe reaction.In this paper we have investigated the thermal decayrates of trajectories very close to the NHIM based onequilibrium trajectories located exactly on the NHIM. Infuture work it will be necessary to also study trajectoriesout of the NHIM. An important question is whether onecan define a thermal equilibrium or at least a stationarydistribution on the DS in these nonequilibrium systemswith which one can obtain the reaction rate.Recently, we have investigated the influence of exter-nal driving on decays in the geometry of the LiCN iso-merization without considering noise and friction [55].Meanwhile, the dissipation arising from an argon bathon that reaction was seen to be representable by theLangevin terms [34]. Thus, the methods presented hereopen up the possibility of considering the thermal effectsin a driven LiCN isomerization reaction and other chem- ical reactions of interest.
V. ACKNOWLEDGMENTS
The German portion of this collaborative work wassupported by Deutsche Forschungsgemeinschaft (DFG)through Grant No. MA1639/14-1. RH’s contribu-tion to this work was supported by the NationalScience Foundation (NSF) through Grant No. CHE-1700749. M.F. is grateful for support from the Landes-graduiertenf¨orderung of the Land Baden-W¨urttemberg.This collaboration has also benefited from support bythe European Union’s Horizon 2020 Research and Inno-vation Program under the Marie Sk lodowska-Curie GrantAgreement No. 734557. [1] H. Eyring, J. Chem. Phys. , 107 (1935).[2] E. P. Wigner, J. Chem. Phys. , 720 (1937).[3] K. S. Pitzer, F. T. Smith, and H. Eyring, The TransitionState , Special Publ. (Chemical Society, London, 1962)p. 53.[4] P. Pechukas, Annu. Rev. Phys. Chem. , 159 (1981).[5] B. C. Garrett and D. G. Truhlar, J. Phys. Chem. ,1052 (1979).[6] D. G. Truhlar, A. D. Issacson, and B. C. Garrett, “The-ory of chemical reaction dynamics,” (CRC Press, BocaRaton, FL, 1985) pp. 65–137.[7] G. A. Natanson, B. C. Garrett, T. N. Truong, T. Joseph,and D. G. Truhlar, J. Chem. Phys. , 7875 (1991).[8] D. G. Truhlar, B. C. Garrett, and S. J. Klippenstein, J.Phys. Chem. , 12771 (1996).[9] D. G. Truhlar and B. C. Garrett, J. Phys. Chem. B ,1069 (2000).[10] T. Komatsuzaki and R. S. Berry, Proc. Natl. Acad. Sci.U.S.A. , 7666 (2001).[11] H. Waalkens, R. Schubert, and S. Wiggins, Nonlinearity , R1 (2008).[12] T. Bartsch, J. M. Moix, R. Hernandez, S. Kawai, andT. Uzer, Adv. Chem. Phys. , 191 (2008).[13] S. Kawai and T. Komatsuzaki, Phys. Rev. Lett. ,048304 (2010).[14] R. Hernandez, T. Bartsch, and T. Uzer, Chem. Phys. , 270 (2010).[15] O. Sharia and G. Henkelman, New J. Phys. , 013023(2016).[16] J. C. Keck, Adv. Chem. Phys. , 85 (1967).[17] R. G. Mullen, J.-E. Shea, and B. Peters, J. Chem. Phys. , 041104 (2014).[18] W. H. Miller, J. Phys. Chem. A , 793 (1998).[19] M. Feldmaier, A. Junginger, J. Main, G. Wunner, andR. Hernandez, Chem. Phys. Lett. , 194 (2017).[20] 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).[21] E. Pollak, J. Chem. Phys. , 1116 (1990).[22] T. Uzer, C. Jaff´e, J. Palaci´an, P. Yanguas, and S. Wig-gins, Nonlinearity , 957 (2002). [23] T. Bartsch, R. Hernandez, and T. Uzer, Phys. Rev. Lett. , 058301 (2005).[24] S. Wiggins, Regul. Chaotic Dyn. , 621 (2016).[25] A. Junginger, G. T. Craven, T. Bartsch, F. Revuelta,F. Borondo, R. M. Benito, and R. Hernandez, Phys.Chem. Chem. Phys. , 30270 (2016).[26] G. T. Craven, A. Junginger, and R. Hernandez, Phys.Rev. E , 022222 (2017).[27] R. Brown, Philos. Mag. , 161 (1828), ibid. , 161 (1829).[28] H. A. Kramers, Physica (Utrecht) , 284 (1940).[29] P. H¨anggi, P. Talkner, and M. Borkovec, Rev. Mod.Phys. , 251 (1990), and references therein.[30] V. I. Mel’nikov and S. V. Meshkov, J. Chem. Phys. ,1018 (1986).[31] E. Pollak, H. Grabert, and P. H¨anggi, J. Chem. Phys. , 4073 (1989).[32] P. L. Garc´ıa-M¨uller, F. Borondo, R. Hernandez, andR. M. Benito, Phys. Rev. Lett. , 178302 (2008).[33] P. L. Garc´ıa-M¨uller, R. Hernandez, R. M. Benito, andF. Borondo, J. Chem. Phys. , 204301 (2012).[34] A. Junginger, P. L. Garc´ıa-M¨uller, F. Borondo, R. M.Benito, and R. Hernandez, J. Chem. Phys. , 024104(2016).[35] A. J. Lichtenberg and M. A. Liebermann, Regular andStochastic Motion (Springer, New York, 1982).[36] R. Hernandez and W. H. Miller, Chem. Phys. Lett. ,129 (1993).[37] R. Hernandez, Ph.D. thesis, University of California,Berkeley, CA (1993).[38] E. Ott,
Chaos in dynamical systems , second edition ed.(Cambridge University Press, Cambridge, 2002).[39] S. Wiggins,
Normally hyperbolic invariant manifolds indynamical systems , Vol. 105 (Springer Science & BusinessMedia, 2013).[40] M. Feldmaier, R. Bardakcioglu, J. Reiff, J. Main, andR. Hernandez, J. Chem. Phys. , 244108 (2019).[41] J. S. Langer, Ann. Phys. (N.Y., NY, U.S.) , 258 (1969).[42] E. Eli Pollak, J. Chem. Phys. , 973 (1995).[43] L. Farkas, Z. Phys. Chem. (Leipzig) , 226 (1927).[44] E. Pollak and P. Talkner, Chaos , 026116 (2005).[45] M. Tsch¨ope, M. Feldmaier, J. Main, and R. Hernandez,Phys. Rev. E , 022219 (2020). [46] T. Bartsch, T. Uzer, and R. Hernandez, J. Chem. Phys. , 204102 (2005).[47] T. Bartsch, F. Revuelta, R. M. Benito, and F. Borondo,J. Chem. Phys. , 224510 (2012).[48] J. Reiff, R. Bardakcioglu, M. Feldmaier, J. Main, andR. Hernandez, “Enhancing and controlling rates in modelchemical reactions through external driving,” in prepa-ration.[49] R. Bardakcioglu, A. Junginger, M. Feldmaier, J. Main,and R. Hernandez, Phys. Rev. E , 032204 (2018).[50] R. Kubo, Rep. Prog. Phys. , 255 (1966).[51] J. Keizer, J. Chem. Phys. , 1679 (1976).[52] R. Hernandez and F. L. Somer, J. Phys. Chem. B ,1064 (1999).[53] P. Schraft, A. Junginger, M. Feldmaier, R. Bardakcioglu,J. Main, G. Wunner, and R. Hernandez, Phys. Rev. E , 042309 (2018).[54] M. Kuchelmeister, J. Reiff, J. Main, and R. Hernandez,Regul. Chaotic Dyn. , 496–507 (2020).[55] M. Feldmaier, J. Reiff, R. M. Benito, F. Borondo,J. Main, and R. Hernandez, J. Chem. Phys. , 084115(2020).[56] G. T. Craven, T. Bartsch, and R. Hernandez, J. Chem.Phys. , 041106 (2014).[57] R. Zwanzig, Phys. Rev. , 983 (1961).[58] See Supplemental Material athttps://link.aps.org/supplemental/10.1103/PhysRevE.102.062204for average decay rates at various temperatures. upplemental Material for “Thermal decay rate of an activated complex in a drivenmodel chemical reaction” Robin Bardakcioglu, Johannes Reiff, Matthias Feldmaier, J¨org Main, and Rigoberto Hernandez ∗ Institut f¨ur Theoretische Physik 1, Universit¨at Stuttgart, 70550 Stuttgart, Germany Department of Chemistry, Johns Hopkins University, Baltimore, Maryland 21218, USA (Dated: January 6, 2021)
This supplemental material includes the stroboscopicmaps of the thermal trajectories on the NHIM for γ = 0 . ω = π (top) and ω = 0 . π (bottom)at various temperatures other than T = 0 . T = 0 .
1, 0 .
3, 0 . .
7, respectively. Asclaimed in the main text, the “moat” is visible in all ofthe lower frequency driving cases. ∗ Correspondence to: [email protected] a r X i v : . [ phy s i c s . c h e m - ph ] J a n − . − . . . . y − − − v y (a) . . . . h k i − . − . . . . y − v y (e) h k i − . − . . . . y − − − v y (b) . . . . h k i − . − . . . . y − v y (f ) h k i − . − . . . . y − − − v y (c) . . . . h k i − . − . . . . y − v y (g) h k i − . − . . . . y − − − v y (d) . . . . h k i − . − . . . . y − v y (h) h k i FIG. 1. Stroboscopic view of a thermal trajectory at (a, e) T = 0 .
1, (b, f) T = 0 .
3, (c, g) T = 0 .
5, and (d, h) T = 0 . γ = 0 . ω = π (a–d) and ω = 0 . π (e–h). The color maps show the time-averaged decayrates h k i of a non-thermal activated complex on the NHIM for the initial conditions ( y, v y ) at time t0