rameter, c4 . conditions and noise sequences from the appropriate ensem‡ e factor Z in Vα trajectory = v(0) = until (λu −its λsen)u(0). (12) In order to calculate the corresponding reaction rate, k, bles, and simulate the behavior of each e motion.Bartsch It et al. 4510-4 J. Chem. Phys. 136, 224510 (2012) the simplest DS, defined by x = 0, and use the basic we choose ergy is so far below the barrier top that it(10) can be regardedonly as a formal solution, In general, Eq. represents 091102-2 et al. J. J. Chem. Phys. 091102-2 Revuelta Revuelta et al. Chem. Phys.136, 136,091102 091102(2012) (2012) flux-over-population rate formula (see e.g., Ref. 20), which having been thermalized on either the reactant or the product since itscoordiright-handwhere side depends on the unknown functions here q is an N dimensional vector of mass-scaled states that k is proportional to the reactive flux side of the barrier. The trajectory can then be classified as re"u and "s. However, since it is known that the critical trajec-1 tes, The U(q) isdynamics the potential of mean force, ! isdefined a symmet1 dynamics of a unit mass particle defined by coordiThe of a unit mass particle by coordivv v v ‡state it reached. All ‡ active or non-reactivetory depending on what (t) = S[λ ; t], a saat(t)all= times, − S[λs , ξα ; t], b b uneighborhood u , ξαbarrier react remains in the of the positive-definite N × N matrix of damping constants, and nate x moving in in anumerical one-dimensional potential U(x) is given nate x moving a one-dimensional potential U(x) is given (17) J = %v&α,v λu − λs in this λu −λs results presented in this work were obtained reactive reactive (t) is we the have fluctuating force exerted by the heat It is con- force will be small, and thebath. anharmonic then its influence (17) by by and way. 30, 41 cted to the friction matrix ! by the fluctuation–dissipation can be evaluated through perturbation theory, thus obtaining are given by and the S functionals across the DS. This flux is to be averaged over different realcular, 40we did non non " algorithm " This is conceptually straight-forward, but com‡ ‡  2 *‡ (x)(x) −− γ x˙γ + ξα ξ(t). x¨ = −U x˙ + (t). (1) x¨ = −U ‡(1) eorem α reactive reactive c4 V1 + c4 V2 ∞+ . . . 2in powers of the an expansion Vα 1= V0 + x x izations of the α, and also over a Boltzmann ensemble x noise, x 3 1  putationally costly. It would be highly,desirable to find a cri-g(τ ) exp(µ(t ! "  Fabio Revuelta Thomas Bartsch , Florentino Borondo , Rosa M. Benito − − τ )) dτ : Re µ > 0,  T " " anharmonic coupling parameter c . of starting at the DS. No4  isξ αthe fluctuating force exerted by the bath, which Here, ξ α ξ(t)α (t) is(t)ξ the fluctuating force exerted by the bath, which Here, 091102-2 Revuelta et al. J. Chem. Phys. 136, 091102 (2012) 091102-2 Revuelta et al. J. Chem. Phys. 136, 091102 (2012) ‡ ‡initial velocities, v, for trajectories (t ) = 2k T ! δ(t − t ), (9) B terion that allows one to identify the reactive trajectories withVV α t α uu In the harmonic approximation, the critical trajectory is 1 S [µ, g; t] = ‡ ‡ ‡ ‡ tice thatPolitécnica only reactive trajectories should be136, included in the *Física s connected to to thethe damping strength, γout ,Sistemas the fluctuation– is connected damping strength, γby , aby the fluctuation– τDepartamento Grupo de Complejos and de y Mecánica. E.T.S.I.A. Universidad de Madrid. 28040 Madrid (Spain) 224510-8 Bartsch etThe al. following J. Chem. Phys. 224510 (2012) x t ,v t x t ,v t t out having to carry numerical simulation.  is the Boltzmann constant and T is the temperature. here k  B 2 School of Mathematics, Loughborough given by  dissipation theorem dissipation theorem average. University, LE11 g(τ ) exp(µ(t −coordiτ )) dτ : Re µ <3TU 0. (United +particle The dynamics of a unit mass defined bybyLoughborough coordiThe dynamics of a unit mass particle defined v v v v Kingdom) sections will describe the phase space structures that will prohroughout most of this work, we3 consider a one-dimensional a aassumption b b can −∞(ICMAT), Universidad Autónoma de The non-recrossing of conventional TST Departamento de Química, and Instituto de Matemáticas Madrid, Cantoblanco, 28049 Madrid (Spain) nate x moving in a one-dimensional potential U(x) is given nate x moving in a one-dimensional potential U(x) is given "such " " " ‡ λs t vide a criterion. Eq. (3) being absent in one dimension—to find s s (18) -0.236 (t)ξ (t )$ = 2k T γ δ(t − t ). (2) #ξ (t)ξ (t )$ = 2k T γ δ(t − t ). (2) #ξ oblem in which the friction matrix ! simply reduces to a "u (t) = 0 and "s (t) = −x (0)e , (13) 0 0 eaction. The α α α α α α B B reactive reactive be restated here by saying that the reactive trajectories are by by %& # $ The subscript τ is used in the S functional to indicate the inte‡2 alar γwell , andunthe position vector q contains a single coordinate n is V The model potential that we have chosen to study is The model potential that we have chosen to study is those crossing the DS with velocity v > 0. We call the rate non non FIG. 1. Schematic view of the phase space structure near the transition state FIG. 1. Schematic view of the phase space structure near the transition state and for this case Eq. (12) yields -0.237 " " , (40) κ = exp − variable. This subscript will be left out whenever this A BSTRACT If we expand the potential of mean force around its saddle ˙ ˙ (x) − γ x + ξ (t). (1) x¨gration −U (x) − γ x + ξ (t). (1) x¨==trajectory −U the transmisα α TST 2k T forfor thethe Langevin equation forfor thethe harmonic (a)reactive and (b)(b) with this trajectory Langevin equation harmonic (a) andanharmonic anharmonic reactive B x x . Any effects constant obtained approximation k x x α III. TIME-DEPENDENT INVARIANT MANIFOLDS does not cause any ambiguities. Similarly, we have for the int, we can write it as 1 c 1 c In this work, we report that allows the reactive trajectories exactly in a system described by an anharmonic potential that 4 44 method ‡ identification ‡ ‡ ofinvariant 2 2 2a new 4 in this work. cases. The time-dependent manifolds are attached to the TS trajeccases. The time-dependent invariant manifolds are attached to the TS trajec20 V(3) ≡ Vof0 simplicity = (λ λnot (0)bybythe U (x) == − −ωbωxb x+ + x x, Here, (3) U (x) ,Here,ξξαα(t) ‡ -0.238 u − exerted s )uindicated beyond TST are customarily summarized a transmisis the fluctuating force bath, which (t) is the fluctuating force exerted the bath, which αsake ‡ ‡ into This where only the average over the noise remains. expres(t) and in our notation that u V V tory, and move through phase space with it. In the harmonic limit, they appear tory, and move through phase space with it. In the harmonic limit, they appear interacts with its environment without the need of any numerical simulation. The method is based on the identification of the geometrical structures (the 2 4 2 4 nnormalized u u c c3 3 1 2 2 TST 4 ‡ damping ‡ ‡ . forxa‡ xt‡harmonic sion coefficient κ =ink/k is,isconnected totothe strength, γγ(a), , ,but byby the fluctuation– connected the damping strength, the fluctuation– as (red) straight lines (a), they get deformed by anharmonic couplings (b). as (red) straight lines but they get deformed by anharmonic couplings (b). s (t) depend on the realization α of the noise, although they sion was derived Ref. 33 barrier. It is now x 4model εThex Langevin + ε2 that + · · · (10) U (x) =invariant − ωb xA.+manifolds) ,v t ,vt t by act as separatrices of phase space. Furthermore, we have succeeded in computing exact reaction rates computing analytical that was already derived in Ref. 16. When the solution (13) is although ourour any onealthough derivation equallyapplies todissipation anyother other one2derivation 3equally 4 appliestodissipation -0.239 and non-reactive trajectories areare represented inin black. Reactive and non-reactive trajectories represented black. theorem theorem In terms of our stochastic invariant manifolds, reactive bothReactive obviously do. clear that the same expression holds also for anharmonic des the entirecase. corrections to the famous Kramers’ formula. substituted into Eq. (10), x = x‡ + "u + "s is replaced by dimensional This generality will be emphasized by usdimensional case. This generality will be emphasized by us‡ ‡ ‡‡ , as discussed before. here ε is a formal perturbation parameter that serves only trajectories are characterized by v > V (t) and s (t) solve the equations of moThe functions u potentials if the critical velocity V α is suitably modified. "" " " 3 3 s s armonic bar(t)ξ (t )$ = 2k T γ δ(t − t ). (2) #ξ (t)ξ (t )$ = 2k T γ δ(t − t ). (2) #ξ We begin by specifying the model that will be used. The -0.240 ng f(x), equal to −c x in our case, to denote any anharmonic ing f(x), equal to −c x in our case, to denote any anharmonic α α B α α B α α 4 4 ‡ harmonic λs t ‡ limit f (x) = 0. They can therefore keep track of the orders of perturbation theory, and finally X(t) tion=in the be criterion, Using this the Boltzmann average overinvelocities in Remarkably, no anharmonic corrections arise the rate exx (t) − e x (0), (14) -10 -5 0 5 10 α α is to evaluate Langevin equation describes the reducedregarded dynamics of a lowforce. The to to higher dimension is is straightforward force. The higher dimension straightforward ll be set to εextension =extension 1. For the mean force itself we write asλhave the coordinates of0,'u a'u special trajectory called the The that we chosen to study is Themodel modelpotential potential that we have chosen to study is c Eq. (17) can be evaluated, as it was in Ref. 16 the harmonic 1. MODEL: LANGEVIN EQUATION pression (40). 3. REACTION RATE CALCULATION FIG. 1. Schematic view ofof thethe phase space structure near thefor transition state FIG. 1. Schematic view phase space structure near the transition state 15 Since > 0 and λ < 0, increases exponentially in time, Since λ > 0 and λ < increases exponentially in time, 4 4 u s u s ved by a nudimensional system coupled to an external heat bath. It is andand will be be presented elsewhere. will presented elsewhere. 224510-4 Bartsch etLangevin al. J. C which is the harmonic approximation to the coordinate x(t) of trajectory for the equation for the harmonic (a) and anharmonic (b) trajectory for the Langevin equation for the harmonic (a) and anharmonic (b) TS trajectory. This trajectory is distinguished from all other dU Ifthe we have a perturbative expansion case. This gives the exact expression whereas 's shrinks accordingly. More importantly, lines whereas 's shrinks accordingly. More importantly, the lines 2 1 c 1 c The transmission factor is calculated by making theare average over the noise [6] o thisForcrucial 4 4 4described given The dynamics of low-dimensional system coupled to an external heat bath can be every fixed realization of the noise, the LE gives rise For every fixed realization of the noise, the LE gives rise =aby ω x + f (x), (11) − 2 2 2 2 4 cases. The time-dependent invariant manifolds attached to the TS trajeccases. The time-dependent invariant manifolds are attached to the TS trajecb trajectories that are exposed to the noise the critical Equation (10a) then gives FIG. 3.trajectory. CriticalU velocity for one realization of the noisesame for (3) a(3) one- by the fact ω x (x) = − x + , ω x U (x) = − x + , dx b b 13,13, 14 14 vector of mass-scaled 'u 0coordi=4stochastic dynamics of the 'u== 0and and's's =0 0are areinvariant invariant under dynamics of thephase ‡ with ‡it. ‡ where q [1-3] istrajectory. an N dimensional driving under force ⇠the (t). This time dependence can be ‡ et 2harmonic *+ ( ) usingtrajectory Langevin equation ↵the tory, and move through space it. In the appear tory, and move through phase space with In the harmonic they appear no 2 4 2 This orbit reThis orbit reo atotrajectory specific called TS a specific trajectory called TS trajectory. 224510-10 Bartsch al. where q is an N dimensional vector of mass-scaled coordiwhere +‡2ε V2 + · ·limit, ·limit, , they (41) V = V0 + εV1 V dimensional barrier with quartic anharmonicity, c4 , for of ωb = 1, saddle γ = 2.5,point for all that it remains in the vicinity the nates, U(q) is q the potential of=mean force, is a1symmetremoved by thestable coordinate shift as straight lines (a), but they getget deformed anharmonic couplings (b). as(red) (red) straight lines (a), but they deformed by anharmonic couplings (b). αby system; being the unstable and manifolds of the system; being the unstable and stable manifolds of theorigin, origin, here f (x) denotes the anharmonic parts of the force. ¨ ˙ uence, every = r U(q) q + ⇠ (t), (8) nates, U(q) is the potential of mean force, ! is a symmetk T 1. Numerical simulation results (red crosses), harmonic approximamains in in thethe vicinity of of thethe energetic for all times, mains vicinity energeticbarrier barrier for all times, q B κ = exp − , (18) 1 ↵ although our derivation equally applies to any other onealthough our derivation equally applies to any other one. times, whereas a typical trajectory would descend into ei‡ "u (t) = S[λ , f (X); t], Reactive and non-reactive trajectories are represented in black. Reactive and non-reactive trajectories are represented in black. ric positive-definite N ⇥ N tion matrix of damping constants, and 1horizontal u As S[λu , ξα ; t], s u (t)to=obwe can substitute into (40) and expand the exponential 2k T respectively. As shown in Fig. 1(a), invariant manifolds respectively. shown in Fig. 1(a), these invariant manifolds (23) (gray line), perturbative results to first-order (23)these + (37) B ric positive-definite N × N matrix of damping constants, and α ‡ ‡ without ever descending into any of the potential wells. Other without ever descending into any of the potential wells. Other λgenerality −Itof λiss conλu − λs dimensional case. This will be emphasized by dimensional case. This will be emphasized usugenerality ther the reactant or the product well both in remote uby=usu the u ,C. s past = s for s , the one-dimensional (16) Results potential (t) mass-scaled is the fluctuating force exerted by the heat bath. 0.37 where q are the⇠↵N coordinates, U(q) is the potential mean force, Γ is a (green straight line) and second-order (23) + (37) + (38) (blue line). tain a series of rate corrections separate trajectories qualitative behavior: Those trajectorieswith withdifferent different qualitative behavior: Those (t) is the fluctuating force exerted by the heat bath. It is conξ αperturbation Using theory, we finally get [7-8] 3 3separate rajectories will be referred to the TS trajectory, which will trajectories will be referred to the TS trajectory, which will ing f(x), equal to −c x in our case, to denote any anharmonic ing f(x), equal to −c x in our case, to denote any anharmonic nected tostates theNxN friction matrix by the fluctuation–dissipation and in the distant future. Accordingly, when using coordi4 4 symmetric positive-definite matrix of damping constants and ξ (t) is the fluctuating force Time-dependent transition α 30, 41 where only the average over the noise remains to be done. By and therefore the leading-order velocity correction is nected to the friction matrix ! by the fluctuation–dissipation above the stable manifold (larger relative velocity) move to above the stable manifold (larger relative velocity) move to are0.37 giv and the S functionals 35 With the help of Eq. (50) the leading term in the trans2 where Z be be taken as a moving coordinate origin. Since the LE is a taken as a moving coordinate origin. Since the LE is a , theorem 8 1 force. The extension to higher dimension is straightforward force. The extension to higher dimension is straightforward κ = κ0 + εκ + · the · · , critical velocity (42)  * nates (utheorem and (s, we are describing a trajectory relative to the perturbative exerted by the heat bath, which satisfies the fluctuation-dissipation 40 1 + ε κ2 for > theorem substituting expansion &ininthe ' ‡ > 3 thetheproduct side distant future, those below itit product side the distant future,while while those below ∞ > Since λ > 0 and λ < 0, 'u increases exponentially in time, Since λ > 0 and λ < 0, 'u increases exponentially in time, mission factor (43a) can be evaluated, giving Because the differential Langevin equation (8) Dis second-order dif‡ g(⌧) exp(µ(t ⌧)) d⌧ : Re µ > 0 E 3 V = c S [ , X ; 0] > u s u s which is again consistent with Eq. (30), and second order equation, itsaitsphase space is twosecond order differential equation, phase space is two4 u > be presented elsewhere. and will be= presented elsewhere. the TS trajectory, which acts as a moving coordinate origin. 2 > 1 1 V S[λ , f (X); 0] = −c S λ , X ; 0 (15) T and 0 will 0 ! expression " factor  u t ).move, 4other u hand, > ‡ ‡ this t famous − g(τ ) exp( 1 where κ is the Kramers’ [9] into and expanding the exponential, a pertur < T ⇠ (t)⇠ (t ) = 2k T (t (9) " " on the other hand, towards the reactant side. move, on the towards the reactant side. where 0.37 B 0 u (t) = S [ , ⇠ ], s (t) = S [ , ⇠ ], ↵ ↵  rential equation, its phase space is two-dimensional, with whereas 's shrinks accordingly. More importantly, the lines whereas 's shrinks accordingly. More importantly, the lines t u ↵ t s ↵ ξ (t)ξ (t ) = 2k T ! δ(t − t ), (9) ˙ x. ˙If If dimensional, with coordinates x and v= x. we now introdimensional, with coordinates x and v= we now intro↵For S [µ, g; t] = B > ⌧ what follows, we will refer (u and (s α fixed ofofthe the LE gives rise Forevery every3c fixed realization thenoise, noise, theto LE gives riseas>>relative α 2In realization Z t coort α 2 λ u s u s u > + c κ + c κ + . . ., bative series of rate corrections, κ = κ ‡ 1 In coordinates, the invariant manifolds appear In space coordinates, the invariant manifolds appear S g; t] = 2space λfixed 31414 3 4 > 0 the 4 1 2[µ, ˙ As it was observedfor * s τfixed intoRefs. 30–32, ordinates x new and vcoordinates τ 13, 13, 'u = 0 and 's = 0 are invariant under the dynamics of the 'u = 0 and 's = 0 are invariant under dynamics of the , 4 > = . (53) κ x = x. duce thethe new coordinates duce the quartic potential (3). Given the initial condition t (17) > 0 = − trajectory Sτ [λ , X aTS (τTS )(e S[λu , Xu and ; 0]−S[λ Xre;vτx] :as+kspace Vspecific fixed dinates to the original s, or This orbit aaspecific called trajectory. This orbit retrajectory called trajectory. g(⌧) exp(µ(t ⌧)) d⌧ : Re µ0 < 0is  u and ux,and $P % κ = , (43a) 4work, Throughout most oftothis we will study oneis the Boltzmann constant and T the temperature. where α  0.36 B ω 26,36  is obtained, where b to the TS trajectory, as shown by the dashed lines in attached to the TS trajectory, as shown by the dashed lines in ‡ap-λu −λs attached 1 e dynamics2.ofTHE the Langevin equation in the harmonic system; being system; beingthe theunstable unstableand andstable stablemanifolds manifoldsofofthe theorigin, origin, g(τ ) exp + and the Sfor functionals are given by 1 "s (0) = x (0) − "u (0), "s (t) can be obtained from Eq. Z coordinates. mains in the vicinity of the energetic barrier all times, mains in the vicinity of the energetic barrier for all times, TIME-DEPENDENT TRANSITION STATE AND THE 1 1 1 problem in which the position vector q contains a x x x x −v + λ −v + λ v− λdimensional v− λ 1 α Throughout most of this work, we consider a one-dimensional s s u u −∞ ⌧ 1 and κ , with i≥1, are corrections due to the anharmonicities, which are zero for odd i. ‡ Fig. 1(b); their instantaneous position depends on the realizaFig. 1(b); their instantaneous position depends on the realiza3 oximation canube= by itIf in coordinates 3 , ,coordinate s= , (4) srewriting = x. (10b), , (4) udiagonalized = single i respectively. As shown in Fig. 1(a), these invariant manifolds respectively. As shown in Fig. 1(a), these invariant manifolds ¯ ZOther This is the famous Kramers result for the transmission 1 V2 = c4 X (⌧) e d⌧ 0.36 we expand the potential of mean force The equations of motion in relative coordinates are 8 and then substituted into Eq. (10a) to find "u (t). In , X ; τ ]); 0]. (38) + S[λ without ever descending into any of the potential wells. Other without ever descending into any of the potential wells. 1 problem in which the friction matrix ! simply reduces to a 2 s Z 1 COMPUTATION OF THE INVARIANT MANIFOLDS λ − λ λ − λ λ − λ λ − λ u u s s u u s s 1 ‡ ‡ > κ = %E& , (19a) 0 17 tion of the noise. Accordingly, the manifolds move through tion of the noise. Accordingly, the manifolds move through ‡ > and s given by ‡ 0 α 1 separate trajectories with different qualitative behavior: Those separate trajectories with different qualitative behavior: Those > around its saddle point, we can write it as κ $P V = − V % , (43b) ‡ g(⌧) exp(µ(t d⌧s )u: (0) Re µ >1q0,contains a single 2 subscript factor. V =the ( ⌧)) The τ is used in the S fu α =c > u⌧ scalar γ , and vector coordinate will totothe TS will trajectories willbebereferred referred thecorrection TStrajectory, trajectory,which which will thistrajectories way, the second-order velocity uposition 0 1 > 0 V X (⌧) e d⌧ 3 f (x) > Z4.and 1 DEGREE 1 kB TONE 8(37) > ON A SYSTEM OF OF FREEDOM 1t RESULTS Z 1 A comparison of phase the perturbative corrections (38) space but they still separate trajectories with differphase space but they still separate trajectories with differwith with < 8 above the stable manifold (larger relative velocity) move to above the stable manifold (larger relative velocity) move to 0 gration variable. This subscript w ˙ (u + , (19a) ( u = λ > x. If we expand the potential of mean force around its saddle The perturbation expansion is set up in such a way that u − λ x + λ x −v v be taken as a moving coordinate origin. Since the LE is a be taken as a moving coordinate origin. Since the LE is a > S [µ, g; t] = (18) > Z c c x s x u ‡ > 3 > ⌧ 3 4 8 > g(⌧) exp(µ(t ⌧)) d⌧ : Re µ > 0 V = c4 S [ u , X ; 0] Z 12 3 equation 2 4results 1 , 2 2 with > > Inu = a system of one degree of freedom, Langevin can be solved using a set of > λ − λ t g(⌧) exp(µ(t ⌧)) d⌧ : Re µ > 0 > > u s numerical is shown in Fig. 3. Again, this com, s = (12) > ent asymptotic behaviors. The stable manifold intersects the ent asymptotic behaviors. The stable manifold intersects the U(x) = ! x + " x + " x + . . . . (10) > > > 3c CRITICAL AND ambiguities the product side ininthe distant future, while those below itnot the product side the distant future, while those below itVELOCITY ! " ! " # $ # $ > does cause any point, we can write it as 1 1 ‡ > bsecond > > ‡ ' ( 32of 2 t Z > 1/2 1/2 2 λ τ 3 4 1 1 effectively the noise carries a factor of ε. The critical velocity < > order differential equation, its phase space is twosecond order differential equation, its phase space is twos > 8 > 2 2 2 2 g(⌧) exp(µ(t ⌧)) d⌧ : Re µ > 0 Let us make a perturbative expansion a nonlinear potential as [7] ‡2 ‡2 ‡ ‡ 1 t 3 = attached 4 to S λλcoordinates −= λs − − −+ λ4ω V2⌧))=d⌧ c:4 SRe [ µ X 0. ; 0] > > < relative the so called transition state uu,s u+that smove randomly, u, < − [λ , X (τ )(e S[λ , X ; 0] V TRANSMISSION FACTOR 2 ‡ ‡ + g(⌧) exp(µ(t > : γ (± γΔu,±γΔλs) γ> , (5) , (5) λu,s = 4ω S [µ, g; t] = τ u u > > ⌧ $P V P V = V − V % κ 2 parison confirms the accuracy of the perturbative results. Z . Trajectories with ini. Trajectories with inix = 0 axis in a point with velocity V x = 0 axis in a point with velocity V ‡ S [µ, g; t] = b b 2 > 2other α c of33 simplicity > sake not indicated ‡ = on move, the hand, towards reactant side. move, on other hand, towards reactant side. > ‡α t− λ > ⌧ 3the 0 ‡,the 1 the 0 2 (QUARTIC tintro> Z < α α > 2 > ‡ λ PERTURBATION) V c S [ , X ; 0] ¥ ¥ 1 > g(⌧) exp(µ(t ⌧)) d⌧ : Re µ > 0 2 2 ˙ ˙ dimensional, with coordinates x and v = x. If we now introdimensional, with coordinates x and v = x. If we now t 0.0 > 3 u  = c k T u s Z 1 V = c S [ , X ; 0] c c 1 > is linear in the noise. If V is one order ε higher, it must V 2(k T ) k T 4 B > trajectory sBartsch This particular trajectory is the only solution of the equations of 1 Phys. > 3 4 4 u 8) [4-6]. V = c S [ , X ; 0] ‡ ‡ S [µ, g; t] = B B 2 > 4 ‡ 224510-7(u , where et al. J. Chem. 136, 224510 (2012) 1 2 2 3 2 4 > 3 > 4 u ⌧ > 0 1 2 2 Z t thatX " is a > formal perturbationThe parameter serves only to > 2(10) > f > s appear (t) depend on α +(x) g(⌧) exp(µ(t ⌧)) d⌧ : Re µ < 0 function introduced in Eq. (29) plays a special role !bthe tial position x = 0 and initial velocity v larger than this critical tial position x = 0 and initial velocity v larger than this critical ω x x x + ε + ε + · · · , U (x) = − : > > κ EV = − V , (19b) > s ( realization u s) In space fixed coordinates, the invariant manifolds appear In space fixed coordinates, the invariant manifolds t > + g(⌧) exp(µ(t ⌧)) d⌧ : Re µ < 0 1 b > < -0.236 Z : ‡ ‡ 3 0 1 > > 8 duce the new coordinates duce the new coordinates αConsequently, motion that remains in the vicinity of the saddle point for all times. 3 0 finally¯ will 3= (˙ s λ (s − . (19b) 1 > 2 3 4 (a) > s 1 relative values are defined by the time-dependent shift relative values are defined by the time-dependent shift V = c S [ , X ; 0] g(⌧) exp(µ(t ⌧)) d⌧ : Re µ > k T cubic. be quadratic in the noise, and V Skeep [µ, g; t] = ‡ ‡ > track of the orders of perturbation theory, and ' ( The subscript ⌧ is used in the S functional to indicate the inteZ > > 4 u 1 + g(⌧) exp(µ(t ⌧)) d⌧ : Re µ < 0 B , X ; τ ] + S[λ , X ; τ ]); 0] (16) −S[λ ⌧> 1 both obviously do. : Z s.t 1 u velocity, sare 2 2the‡ un- with > in the perturbation expansion because it represents ‡2 , are reactive, while trajectories with initial veloci, reactive, while trajectories initial velocivelocity, V V λ − λ 3 attached to the TS trajectory, as shown by lines attached to the TS trajectory, as shown bythe thedashed dashed linesin(43c) in > u s -0.237 x = u + s, be> v = λ u + λ (13) α α 0.0 x u s > ‡ > P V − 3 ⌧ 1 V = c S [ , X ; 0] t u ‡ ‡ g(⌧) exp(µ(t ⌧)) d⌧ : Re µ > 0 > set to " = 1. For the mean force itself we write x x −v + λ v − λ x x −v + λ v − λ > < gration variable. This subscript will be left out whenever this 4 u Z where c are the anharmonic terms of the potential. 1 where ε is a formal perturbation parameter that serves only ‡ ‡ ‡ ‡ V = c X (⌧) e d⌧ > s u s u α Z 1 (t) and s (t The functions u 2 1 replaced by 4 ‡ ‡ ‡> ‡ j > 1 2 Z perturbed trajectory. To obtain a different perspective of this 3 are not. The (random) value V therefore encodes are not. The (random) value V therefore encodes ties v < V ties v < V 2k T + g(⌧) exp(µ(t ⌧)) d⌧ : Re µ < 0 > : Fig. 1(b); their instantaneous position depends on the realizaFig. 1(b); their instantaneous position depends on the realizaB S ⌧'u [µ, g; 1 (t), 's = s − s (t), (6) 's = s − s (t), (6) ut] − 'u == u− u=αu(t), ,the (4) uu==At first sight, = , cause (4) -0.238 ,t , ,, , sits"u α 2 (t) = α α=appears α"x 1 ‡1 + > 2 ⌧have ufor that time-dependent and stochas− S[λ , 2c X c X ; t] 0 λutheory, Z Z α α α ‡ ‡ ‡ 3 ⌧ u 3 4 > < does not any ambiguities. Similarly, we the 1 to keep track of the orders of perturbation and finally V = c X (⌧) e d⌧ ‡ 0 t 3 > 3V = c4 tion in the harmonic limit f (x) u⌧ 1 dU X (⌧) e d⌧ 0 1 0 1 he constants λ − λ λ − λ λ − λ λ − λ 1 4 1 κ = − "V V # = > u s u s u s u s λ − λ 2 V = c X (⌧) e d⌧ 8 6 2 2 1 0 function, note that the critical velocity should depend only on ‡ λ t ‡ the relevant information about the realization of the noise conthe relevant information about the realization of the noise conu s S [µ, g; t] = ‡ 4 -0.239 ‡anharmonic s tion of the noise. Accordingly, the manifolds move through tion of the noise. Accordingly, the manifolds move through 2 0 1 with⌧ u (t ⌧) When the Z > 2 0 105 + 830 + 1  system has only a quartic term, c , the B C B C -0.0 Z tic shift (16) has removed both the time-dependence and 1 3 c k T with the abbreviation u (t) = ⇠(⌧) e d⌧ X(t) = x (t) − e x (0). (29) > sake of simplicity not indicated in our notation that u (t) and 0 B C B C 4 B = ! x + f (x), (11) 4 ‡ 0 0 0 c k T ω will be set to ε = 1. For the mean force itself we write ‡ ↵ t $ > 0 B C B C regarded as theCAcoordinates of a B b +the g(⌧) exp(µ(t ⌧)) µ < 0 b d⌧ : Re : B C B where coordinates trajectory are where thethe coordinates of the TSTS trajectory are >  =  V = ( )u (0) ‡ 3 ⌧ @ A @ % 0 u s u s dx u ‡ t 1 #of * ) 2 > (t) for t ≥ 0, but not the behavior of the stochastic force ξ 4 cisely. Inother other words, once the instantaneous of the cisely. words, once the instantaneous position of the 0 4 2 4 phase space but they still separate trajectories with differphase space but they still separate trajectories with differwith with lowest order terms of the critical velocity are given by α Redistribution Downloaded 31 May In 2012 to 141.20.49.4. subject to=AIP license orposition copyright; see http://jcp.aip.org/about/rights_and_permissions -0.240% V c X (⌧) e d⌧ > s (t) depend on the realization ↵ of the noise, although they $ 2 32 dependence on the realization α of the noise. However, this 2c c PHASE SPACE VIEW OF THE TIME-DEPENDENT 4 2 ! 1 +  (1  )( Z ‡2 4 1 > 2 TS trajectory. This trajectory is 2 ‡2 3 3 γ ± (14) = − γ + 4ω λ 0 0 + g(⌧) exp(µ(t ⌧)) d⌧ : Re µ < 0 b -10 -5 0 5 10 Z % % dU 1 . : s,u function represents V − λ ) u (0) (λ 2 1 h of the noise ∞ ∞ the bharmonic This approximation to the S[λ , X ; t] − = − INVARIANT MANIFOLDS OF THE LANGEVIN EQUATION 0ent is a third-order moment and must vanish. Simii u s u c ‡ on the driving at earlier times: Once the initial conditions of a stable manifold is known, any trajectory can unambiguously stable manifold is known, any trajectory can unambiguously ‡ 0 3 0.374 both obviously do. ent asymptotic behaviors. The stable manifold intersects the asymptotic behaviors. The stable manifold intersects the 2 = ω x + f (x), (11) − 1 1 -0.0 ‡ 2 ! " # $ ! " # $ 3 ⌧ ‡ where f (x) denotes the anharmonic parts of the force. is only true in the harmonic approximation. If we express the ‡ ‡ ⌧ ‡ ⌧ 1 1 ‡ u Z b s u ‡ V = ( )u (0) 1 λ (t−τ ) ‡ ‡ λ (t−τ ) . (44) = exp − P = exp − trajectories that are exposed to 1/2 1/2 λ − λ (λ − λ ) (b) u u u s 22 2Z2 t u s V V= u(0) s = 0 d⌧ V = c x (⌧) e x (0) e d⌧ 2 V 1 c X (⌧) e 0 = ( )u 4 dx , , 1 ‡ ‡ 4 (t) = − ξ (τ ) e dτ, (t) = − ξ (τ ) e dτ, u uαcoordinate u s ‡ ‡ 2 , 1 x(t) of the trajectory under study. Moreover, it conα α 0 its −classified γclassified γγ future , (5) λatλu,su,s ==0 be + 4ω − γ‡ ±± its , (5) +1xfrom 4ω 2u all odd-order to Vthe transmission factor must Dv α The functions (t) and scondition (t) solve the equations mo0.372 c4 B T3 2kcorrections T2with 2k trajectory t = are given, fate can only depend on0 xlarly, be the values of initial condition as refrom the values of its initial as re. Trajectories with ini= 0 axis in a point . Trajectories with inix = 0 axis in a point velocity V b b Bwith 0 of velocity ‡ (t ⌧) The transition state trajectory is given by ! 3 ⌧ position coordinate in terms of the relative coordinates (u s HaL α α that it remains in the vicinity e the eigenvalues thatλarise in the diagonalization. They sat = c k T u λ − λ − λ 4 B Z u u s st t s (t) = ⇠(⌧) e d⌧ 2 2 λ τ 2 2 ‡ 2 4 2 3 V = c X (⌧) e d⌧ s ↵ 2 c 1 3They 4 0. 4u ,tial where fbe (x) denotes the anharmonic parts of thevforce. stitutes a suitable basis of the perturbative expansion. tion in the harmonic limit = can therefore be !times, λ × S , X(τ )(e S[λ X ; 0] − S[λ , X ; τ ] 0.370 V1and = The c S [ , X ; 0] 2 s ( u whereas s) zero. According to the fluctuation-dissipation theorem (9), τ u u 3 u  = c k T the future noise. separatrix between reactive and nonreactive or non-reactive. Finally, the presence of anharmonicactive or non-reactive. Finally, the presence of anharmonicreactive u s 1 position x = 0 and initial velocity larger than this critical tial position x = 0 and initial velocity v larger than this critical 1 -0.0 4 B b (s as a typical traject 2 4 3 y λs < 0 < λu and u‡ (t) = (t ⌧) 2 √ u 0 Z ! ( ) and the lowest order corrections to Kramers’ factor become We will now address the problem of evaluating the noise avs u s relative values are defined by the time-dependent shift relative values are defined by the time-dependent shift ⇠(⌧) e d⌧ regarded as the coordinates of a special trajectory called the B. Time-dependent transition states 1 ‡ ‡ The leading-order correction to the critical velocity can 0.368 " b ↵ nonthe noise carries a factor k T , so that a perturbative expan% %t u t 1 s active(ttrajectories must therefore also be determined by only ities(f(x) (f(x)&=&=0)0) will distort invariant manifolds, indiities will distort the manifolds, as indi, are reactive, while trajectories with initial velocivelocity, V , are reactive, while trajectories with initial velocivelocity, V 2theinvariant 2as ther the reactant or the product B Du ¯ α α ‡ ‡1 1 ⌧) + S[λ , X ; τ ]); t . reactive ‡ t u 2 2 (a) s TSxtrajectory. This iserages distinguished from all other Dx (u c4 trajectory in Eq. (43). 3 0.366 + + (s, (20) x = 2 0 ‡ ‡ 1 0 from (24) as ‡ u (t) = ⇠(⌧) e d⌧ ‡ Z ‡uobtained λ (t−τ ) ‡ be λ (t−τ ) + λ = −γ , λ λ = −ω . λ V = c S [ , X ; 0] , 1powers ‡ ‡ ‡ ‡ s s s u s 8 4 in 8the 3 u  = c k T ↵ and distant future. the future noise. Yet the perturbation term in (21) depends, via 2 1 cated by the red lines with arrows in Fig. 1(b). The main step cated by the red lines with arrows in Fig. 1(b). The main step b sion in powers of ε corresponds to an expansion in of are not. The (random) value V therefore encodes ties v < V are not. The (random) value V therefore encodes ties v < V 4 B 6 1 1.5 4 Acco 2 1 c ξ (τ ) e dτ. (7) sα (t) == −− ξ (τ ) e dτ. (7) sα (t) ‡ 3 B. Time-dependent transition states ‡ (t), 's = s − s (t), (6) 'u = u − u (t), 's = s − s (t), (6) 'u = u − u 4 2 4 Z ‡ α α 3 α α α α Z 3 -2 -1.5 -1 -0.5 0 0.5 2 105 + 830 + 1648 + 770 2 1  B C B C trajectories that are exposed to the same noise by the fact that √ 3 c k T 1 α α α α >  = c k T V = ( )u (0) B C B C 4 B 1 u s 4 B 0 ! ( ) 0 0 0 c t u s > V = c S [ , X ; 0] 4 Because the Langevin expression (8) is a second order difh i B C B C s 4u 4 ‡ exp(µ(t ‡ ‡ ⌧)) d⌧‡ : 2 Re c 00 λu λ−u − λs λs −∞−∞ > 3 usthis 2 about 2 B C B C b 3 ‡ theory ‡  =  ⌧ nates (u and (s, we are 2descr g(⌧) µ > 2 u ! ( ) @ A @ A 0 (t), on s (t), which is given by past noise. x > of the to be developed here is the calculation of of the theory to be developed here is the calculation of this the relevant information the realization of the noise conthe relevant information about the realization of the noise con1 k T . By contrast, Eq. (42) is an expansion of the transmiss u s ‡ Ds From (34) we have that 4 ‡ ⌧ ‡ ⌧ = u + s , Eq. (19) turns into with x 4 2 4 4 B V = c X (⌧) e d⌧ > s u b it remains in the vicinity of the saddle point for all times, 3 FIG. 4.10Transm the new set of coordinates, equations of motion read 32 4 !b dif0 e the 12Langevin 0 14equation 1 +  (1  )(3 + + 3) Z 1 Z V = c x (⌧) x (0) e d⌧ ferential equation, its phase space is two-dimensional, with Because (8) is a second-order (a) Critical velocity for a > 8 6 2 2 = > S[λ , f (X); t], "uthe 4 t 0 0 0 0 1 t lead (t) u where the coordinates of the TS trajectory are where the coordinates of the TS trajectory are < the TS trajectory, which acts as 2 0 h i 105 + 830 + 1648 + 770 + 87 1  B C B C 3 c k T If we split up the integration range of the S functional, we deformation. deformation. cisely. In other words, once the instantaneous position of the cisely. In other words, once the instantaneous position of the 3 B C B C 4 B ‡ 1 0 0 T because it has only even-order sion factor in powers of k 2 0 0 0 c whereas a typical trajectory would descend into either the reparticular noise ensemble as a The corresponding equations of motion are The corresponding equations of motion are anharmonicity, c 4 λ − λ B C B C B. Distorted correlation functions B S [µ, g; t] = ‡ 0 u> s ‡ coordinates (t ⌧) ‡ ⌧ ‡ ⌧ ‡ x and v = x ˙ . As it was observed in Refs. 26– B C B C V = c S [ , X ; 0] ⌧ s . s u ferential equation, its phase space is two-dimensional, with  =  ‡ Z Z x 3 x 0(⌧) u@2 A @x (0)2 A e (35)d⌧ 4 % % + (u (s) f (x 1 4+ Z sf↵ (t) ⇠(⌧) e d⌧ 1 V = c e (x) = 1 > t refer toa 4 4 2 > fIn u n what c t i o n follows, o f a q uwe a r t iwill c ∞ ∞ V = − c S[λ , X ; 0], 4 t 32 > strength c for 3 u ! 1 +  (1  )(3 + 10 + 3) actant or the product well both in the remote past and in the This critical velocity can be calculated from the condition This critical velocity can be calculated from the condition find that for t ≥ 0 stable manifold is known, any trajectory can unambiguously stable manifold is known, any trajectory can unambiguously 2 4 1 > terms. 11 harmonic ‡ λvu (u (u˙V= + , (21a) 2 0 As it was observed 0 0 in Refs. 0 30–32, bZ 2 the+dynamics of the Langevin equation in the u˙ = λu uwhich + ‡ 28, ξ (t), u⌧ ˙ = x. coordinates x and v ‡1V ‡ s> ‡uthat anharmonicity ‡ λ (t−τ ) ‡ λ (t−τ ) α x anharmonicity, c , for ω =3.5, γ=2, koriginal T=1. u u 1 from it follows = c S [ , X ; 0] > 0 (t ⌧) 0 1 0 1 = c X (⌧) e d⌧ 3 u dinates and to the u and 4 + g(⌧) exp(µ(t ⌧)) d⌧ : Re µ < 0 s 1 + 'u + 's) f (x + 'u + 's) f (x 3 2 : ‡ ‡ 1 The corrections to the critical velocity that appear in the (t) = − ξ (τ ) e dτ, u (t) = − ξ (τ ) e dτ, u 8 6 4 2 2 λ − λ 0 1 0 1 its Kramers app ! 1 distant future. Accordingly, when using coordinates u and λ − λ λ − λ α α 4 α α u s s (t) = ⇠(⌧) e d⌧ 2 α α ‡ u s u s 105 + 830 + 1648 + 770 + 87 1  be classified from the values of its initial condition as rebe classified from the values of its initial condition as reis contained in is contained in that the trajectory with x(0) = 0 and v(0) = V that the trajectory with x(0) = 0 and v(0) = V 8 6 4 2 2 HbL B C B C (b) Transmission factor for ω =3, γ =7, k T=1. t 3 c k T approximation can be diagonalized by rewritting it in coordi↵ The simplest rate correction can therefore be obtained ‡ B C B C 4 B 3Langevin ⌧ +equation 0 0 the dynamics of the the harmonic αα C (8a) ++ ‡ (8a) u˙ = λu 'u 'u˙'= λu 'u 0in 0 + 87 0 apc105 0 , , Vu− = 0‡ s t tin agreement uFACTOR + 830 1648 + 770 λ λ λ − λ 1 CRITICAL VELOCITY AND (CUBIC PERTURBATION) 4 TRANSMISSION B C B B C B C 1 3 c k T u s B C B C coordinates. B C B C 1 4 B V = c X (⌧) e d⌧ with Eq. (30), and  =  0 0 0 0 0 c simulation result ‡ λ t λ (t−τ ) Numerical simulation results (red crosses), A @ A 0 @ TS 4 relative B C B C u s a trajectory to the trajectory, s s 4s, we are describing averages (43) are expressed in terms of the function X(t), 1 4 λ − λ λ − λ 2 4 2 4 4 2 B C B C  =  u s u s nates u and s given by =the e(30) sstable (0) manifold +manifold e 4ofofthe ξthe ) dτ. In the relative coordinates, oneu , can identify thes (t) reactive f (X); 0]. Vlead = S[λ active ortrajectory, non-reactive. the presence active or non-reactive. Finally, the presence ofanharmonicanharmonicthe find trajectory, stable TSTS this 32 perturbation 0.415 Finally, @trajectory. A(a)To @Tofind Athis 0trajectory. α (τ ! 1 + in (1 of )(3 + (Kramers) 10 +3The 3)= proximation can be diagonalized by rewriting it in coordinates 0, from a quartic the potential. We set c -1.8 4 2 4 4 2 harmoni approximation (gray 0 (b) 0 0 0 0 b horizontal line), equations of motion in r ‡ 32 ! 1 +  (1  )(3 + 10 + 3) which acts as a moving coordinate origin. In what follows, 2 f (x) 1 ‡ 0 (t) and which is in turn given in terms of the components u Z 0 (f(x) 0 the 0 invariantmanifolds, b (s) trajectories uniquely even doing (15) any numerical ‡ 00.410 horizontal as line), first-order correction (green 2and + (u + f (x % % 1 critical ities = & 0) will distort invariant indiities (f(x) = & 0) will distort the manifolds, as indix which will be called critical trajectory, we formally solve the which will be called trajectory, we formally solve the u s given by c ‡without 3 ‡ 4 t t s˙ = λsTo s −obtain − ξ (t). (green line), and 2c which makes V = 0, and calculate the rate correction that α s x)u (0) -2.0and v xucorrections v + x ‡ ‡ we will refer to u s as relative coordinates and to the ‡ 2 ⌧ 3 3 λ τ 2 higher order to the critical velocity line), second-order correction (blue line). x u  = c k T (˙ s = λ (s − . (21b) 1 1 u V = ( ¥ s 1 trajectory. s−cλX 4for Bτ [λu ,3sX(τ sof,‡ a‡ The ‡ in this (t) of the TS They are Gaussian random variables V = (⌧) e d⌧ ‡ (t−τ ) ‡ λ (t−τ ) Vc simulation by the value critical velocity, V , λu calculating − λsf (x λ λ V = S[λ , X ; 0]− S )(e S[λ , X ; 0] 2 4 ‡− ‡u s s 3 0 integral expression depends only on noise = s = , (12) u s 4 u u 2 2 0.405 α 1 cated by the red lines with arrows in Fig. 1(b). The main step cated by the red lines with arrows in Fig. 1(b). The main step equations of motion (8) by equations of motion (8) by V = 0 2 V = ( )u (0) λ − λ Du (t) = − ξ (τ ) e dτ. (7) s (t) = − ξ (τ ) e dτ. (7) s + 'u + 's) f (x + 'u + 's) ˙ (u + ( u = λ u s ! ( ) u s α α u original u and s, or x and v as space fixed coordinates. − λ x + λ x −v v 1 α α 0 s u s . Substituting Eqs. (23) and (37) into (43b), it is is linear in c α we introduce α 0 λ −λ x x s x u in a systematic manner, the expansions c 4 -2.2 u s ‡ Hx HtL,v HtLL b u s u λλs−−λλ 3 ‡, 4 ‡ using perturbation theory [4-5]. Then 2, 0.400 '˙ s = λ '˙ s = λ 's − (8b) 's − , (8b) whose correlation functions were evaluated in Ref. 31. In the u s u s s = , (12) u = −∞ −∞ s s 1 1 t ≥ 0. The term including s (0) contains all the dependence  = c k T V = ( )u (0) of the theory to be developed here is the calculation of this of the theory to be developed here is the calculation of this 4 B 2 2 2 2 2 2 uapproximation, s hese equations decouple in the harmonic i.e., Same as the figure for a quartic The equations of relative coordinates are ‡ 2 motion 3 a 3time42c4 in 0 λuλ− − λ λ U(x, y) = ! x + ! y + cx y λ − λ λ − λ 2 found that u s s 2 u s u s (t) of the TS trajectory represents The position x x y = c k T ¯-2.4 ! ) ‡ ‡ 4 bB s ( notation current and with ‡ 2 ‡on the past, and it drops out when − anharmonicity, c , for a cubic one, c , u s , X ; τ ] + S[λ , X ; τ ]); 0]. (36) S[λ ‡ ‡ 2 2 2 4 0.395 3 or we form X(t). The variable u s 2 deformation. deformation. Ds V = V0 + to ε Vthe + ε V + · · · The corresponding equations of motion are The corresponding equations of motion are f (x) =‡0,‡ but they are still subject time-dependent , f (x + 'u + 's); t] S[λ , f (x + 'u + 's); t] S[λ ) u u -0.3 -0.2 -0.1 0 0.1 0.2 0.3 and ω =1, γ=2, k T=1. s -0.2 -0.1of motion. 0 0.1 !b 0.2s ( u λ t λ t or ‡ u u dependent stochastic driving in these equations ‡ ‡ ‡ ‡ ‡ ‡ 1 ‡ 2‡ ‡ ‡ ‡ in c , which c vu (λ − λ ) c e + , (10a) 'u(t) = C e + , (10a) 'u(t) = C c c f (x) V = 0 (t) = u (t) + s (t) [and v (t) = λ u (t) + λ s (t)]. with x (t) = u (t) + s (t) [and v (t) = λ u (t) + λ s (t)]. with x 4 4 u s u in which the dependence X is the simplest modification of x T γ k u s u s ‡ 3 This critical velocity can be calculated from the condition This critical velocity can be calculated from the condition 4 1 ‡calculated α α velocity α This αforce α α dependence α 2α can the α αgeometrical Bu (0) X (τ )# ; 0]. ochasticα driving ξαα (t). time be (˙s = λs (s − 2 HcL 2 u ˙ = u + , (19a) ‡ S = [λ , "P (54) κ The critical is by computing λ − λ λ − λ ‡ ‡ 2 u τ u α u s u s V = c S [ , X ; 0] Nevertheless, the coordinate shift has removed the stochastic 0 1 0 1 = , (45) σ Not surprisingly, the corrections (35) and (36) depend, "uspace =space ε1 "u "u +vicinity ·vicinity · 3·the 3+ε s,in uthe 4 2ux =8su + ks, T 6v = λ u +4λ s. ++'u ff(x(x 'u++'s) 's) 1+ 2= V c S [ , X ; 0] ‡ (13) ‡ 2 It reads x = u v = u + s. (13) 2== u 2 Now, the geometric phase structure in the of the Now, the geometric phase structure of the x u s α α 2 on past has been removed. 1 x u s is contained in that the trajectory with x(0) = 0 and v(0) V is contained in that the trajectory with x(0) = 0 and v(0) V B moved bystructures the coordinate shift |λ |(λ − λ ) 105 + 830 + 1648 + 770 + 87 1  B C B C ˙ ˙ (the invariant‡ manifolds) that act as separatrices for αα 'u + , (8a) ' u = λ 'u + , (8a) ' u = λ 3 c k T s u s u u B C B C 4 B 0 0 0 0 0 c driving from the leading-order terms in Eq. (21) and pushed through the function X, on the realization of the noise. This 4 B C B C 2 0 B@ 12CAstable 0stable 1A ‡ First, 2harmonic fON (x)manifold 4 SYSTEM B C λ − λ λ − λ barrier can be easily discussed. in the limit, barrier can be easily discussed. First, in the harmonic limit,  =  8 of 6trajectory. 4To 2 this AtOF firstFREEDOM sight, it appears that the 2 u s u s @ A the of the TS find this trajectory, the manifold the TS trajectory. To find trajectory, 0 V = c S [ , X ; 0] 5. RESULTS OF TWO DEGREES "s = − x + ε "s + ε "s + · · · . 0 1 0 1 reactivity. TheseThe structures are linear in the harmonic case and 1 2 3 u 4 4 105 + 830 + 1648 + 770 + 87 s ˙ = s . (19b) 1  4 2 4 4 2 B C B C ‡ constants ‡ 2 3 factBBThe c41k+ constants sT 1 =s−s , 8 potential 6as +010 4 3) 2 inside The average over the noise can be brought the S func2anharmonic C B C B1 it into the anharmonic perturbation. they read, for t ≥ 0, 0 0 32 0 0 dependence reflects the that on an xc4 ! (1  )(3 + C B C , (s (16) (u = u − u 105 + 830 + 1648 + 770 + 87  B C B C c4 both t B C B C tic shift (16) has removed 3 c k T uwill sbe trivially solved, .e.,i.e., f(x)f(x) =g e= 0,t 0, the equations motion areare trivially solved, 0 0 0 0 b  =  B C B C 4 B which called critical trajectory, we formally solve the which will be called critical trajectory, we formally solve the @ A @ A 0 0 0 0 0 0 $ c d ethe f o equations r m e d iof n of t hmotion e (8)(8) 4 4 κ B C B C 4 shape 2 invariant 4 shorthand 4 % 2 q ✓ ◆ # B C B C tional because the latter is notation for an integral. 4 V 1 Z Z 32  =  not only the position, but also the of the manV. CORRECTIONS TO THE REACTION RATES ! 1 +  (1  )(3 + 10 + 3) ‡ ‡ We will write Du 0 @ A @ A 1 ‡ ‡ 1 2 1 2 λ t dependence on the realization α 2 = − s 4 0 of 0study 0systems 0eof more bapply ‡ HtLL ‡ , f (x + 'u + 's); t] S[λ , f (x + 'u + 's); t] S[λ giving presence of anhamonicities. giving 4 2 4 4 2 2 2 One can the developed method to than one degree of freedom. γ ± (14) = − γ + 4ω λ s s equations of motion (8) by equations motion (8) by $s (t)s (0)% = σ , (46a) Hx HtL,v s,u λ t λ t 32 + 'u + 's) f (x + 'u + 's) f (x α ± 2 + 4!ub⌧ (14) s sfirst ! 1 +  (1  )(3 + 10 + 3) b At sight, it appears that the time-dependent and stochas‡ s,u‡2= ⌧ α α κ u 's(t) = C e − , (10b) 's(t) = C e − , (10b) ifolds, are stochastically time dependent. The remaining moment can be evaluated as 0 0 0 0 2 b 0 appr s s V = c X (⌧) e d⌧ is only true in the harmonic 2 V = c X (⌧) e d⌧ '˙ s = λ 's − , (8b) '˙ s = λ 's − , (8b) 3 3 Here, we will show the results obtained for the potential [8] s s 1 for A.λks tGeneral rate expressions λ − λ λ − λ t k +Z λs t 1. λu= t λu1"u "x "s ≥ (31) tic shift (16) has removed both the time-dependence and the u s velocity u s k k λ − λ λ − λ 0 1 , 's(t) = 's(0) e . (9) , 's(t) = 's(0) e . (9) 'u(t) = 'u(0) e 'u(t) = 'u(0) e u We s s Ds calculated the are u also critical numerically for 0 position coordinate x in terms of the eigenvalues that arise in the diagonalization. They satDownloaded 19 Jun 2012 to 158.125.33.105. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/about/rights_and_permissions ‡ ‡ λ I n s y s t e m s are o fthe ore ‡ meigenvalues dependence of the realization ↵ of the noise. However, this that 2 arise in u, f u ⌧the diagonalization. They sat, f (x(x 'u 's); t]t] and (s as S[λ 'u 's); S[λ 1 1 ‡ 3 ‡ ++ 3 ++ 1 1 u u λ t λ t 2 2 2 2 2 2 V = c X (⌧) e d⌧ a given realization of the noise. To this end, an ensemble of 2 2 2 2 2 2 u u ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡ ‡ In a== one-dimensional model, the characteristic function isfy ! λy)s x < 0+<'u(t) λ'u(t) (0) X (τ )# = "u (0) X (τ )# "P u+ Expand the isfy anharmonic as × u xand α 0t e + , ,VELOCITY (10a) = C e + (10a) = C λ U(x, y) = ! y + cx y (t) u (t) + s (t) [and v (t) = λ u (t) + λ s (t)]. with x (t) u (t) + s (t) [and v (t) = λ u (t) + λ s (t)]. with x U(x, = ! ! y + cx y u u s is only true in the harmonic approximation. If we express the dimensions, the 1sstable < 0 3< term and u s u s CRITICAL ‡ ‡ 2 −λ x y x y u α α α α α α α α α α α α u ωαb = λ−uλ− 0 (a) Harmonic approximation of the on the DS λσsλes $u (t)u (0)% , (46b) 2 2 trajectories starting was propagated numerically. 2 2 u− ‡ can be expressed in terms of the critical velocity as χ 2 r position coordinate xthe in terms of the relative coordinates u manifold is not a curve but Now, the phase space structure in the vicinity of the Now, the geometric phase space structure in the vicinity of 2 2 ‡ geometric +( x = x invariant manifolds in relative coordinates. λ 1 λu1+ λ2s 2= −γ , 2 λ u λs = −ωb . 0 u f (X + ε "x1 + ε "x2 ‡++ · · ·)== ε, V f11 + ε0 f2=+ ! · ·2·. ,(b) The " 2 2 2 = 0.4 By recording which trajectories were reactive and which were As is it shown in the figures of this section, (a) λ stable and unstable invariant U(x, y) = ! x + ! y + cx y and s as u s u s V = 0 u Downloaded 31 May 2012 to 141.20.49.4. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/about/rights_and_permissions Downloaded 31 May 2012 to 141.20.49.4. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/about/rights_and_permissions a hypersurface. b discussed. ‡ 2 barrier First, in the harmonic limit, barriercan canbe beeasily easily discussed. First, in the harmonic limit, x y ‡ 1 ‡ (55)‡ ‡ = 3 "u (0) X(τ )# "X (τ )# . (32) are1not, manifolds attached transition : thevtovalue >V , 0 0 2 2 = u + s , Eq. with x x the 0.2 (19) turns the In value ofcould theofbe critical velocity and the of -2 of the critical velocity bracketed with the new set coordinates, the equations motion read ω state trajectory in space-fixed coordinates. i.e., f(x) = 0, the equations of motion (8) are trivially solved, i.e., f(x) = 0, the equations of motion (8) are trivially solved, (39) (v ) = χ b r x ‡ ‡ ‡ ‡ the newr.h.s. set of coordinates, the equations of motion read STABLE MANIFOLD IN Athe SYSTEM . Since f is aswhere termsInon depend on the "x ‡For one fixed j transmission factor approximate to the exact x = x + u + s, (20) high accuracy. realization and for a potential 0 (c) The invariant manifolds get deformed $u (t)s (0)% (46c) V = 0 ‡ ‡α-4 = 0, 0 : v < V . f (x x OF TWO DEGREES OF FREEDOM f (x) 1 , f (x + 'u + 's); t] S[λ , f (x + 'u + 's); t] S[λ giving 1 giving s s due to anharmonicities. λs λ t s t +order. sumed to have an overall order ε or higher, the calculation In the numer ones asterm, weThe increase the perturbation ˙ ( u = λ (u + modified correlation functions that are required here are ˙ u = λ u + ξ (t), u with only a cubic anharmonic the perturbative expanu α 's(t) = C e − , (10b) 's(t) = C e − , (10b) -0.2 s− s λ f (x) 1 ‡ ‡ ‡ λ λu − λs λuλ-6− u s with x = u + s , Eq. (19) turns into λ − λ λ t λ t λ t λ t knowledge of "x for j < k. of the term fk requires only In contrast to the TST approximation (5), and in spite of its u˙ = the u + + ⇠ (t), s u sbe rewritten as a sum u s u s small, but Fig j u ↵ given in Eq. (52). Equation (55) can thus sion is compared to numerical results in Fig. 2. There is good , 's(t) = 's(0) e . (9) 'u(t) = 'u(0) e , 's(t) = 's(0) e . (9) 'u(t) = 'u(0) e CRITICAL VELOCITY AND TRANSMISSION FACTOR -0.4 u s u s -8 Equations (24), (26), and (31) then yield the recurrence relasimplicity, the expression (39) is exact. It allows to evaluate 2λ turbative resu s of exponentially decaying terms, for which the S functional in agreement between perturbative and numerical results. Simi6. CONCLUSIONS ‡ ‡ ‡ 2 −λ t λ t u s f (x) 1 ‡ f (x) 1 -2.4 f (x ˙ + u + $ss)(t)u (a)(0)%α = σ -1 (e -0.5 0− e 0.5 ).1 1.5(46d) 2 2.5 3 f (x 0.04 0 s = λof − noise, − (21a) ξ (t). (15) s˙ = structures ⇠↵ (t). (15) thethat average initial conditions in space Eq. pu⊥ uin s sthe α u˙ =other +Eq. (54) , to evaluate. ss (b) than the firstlar figures are(7)—the obtained for realizations thus (˙ s = λ (s − Thetions calculation of the geometrical act asover separatrices of the phase can be factor is easy This procedure yields y λ + λ (a) s s u λ − λ λ − λ u s u s u s u s u s 0.02 -2.8Obviously, the size of the first -2 used to identify the reactive the need of any numerical simulation. 1 trajectories uniquely withoutDownloaded leading to the same conclusion. For a gen ‡or 31 May 2012 to 141.20.49.4. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/about/rights_and_permissions Downloaded 31 May 2012 to 141.20.49.4. Redistribution subject to AIP license copyright; see http://jcp.aip.org/about/rights_and_permissions 4 2 2 f (x + udecouple + 3c s) 4 σ These equations in the harmonic approximation, i.e., γ The position x‡ (t) of (t) = by S[λ ; t], "ukdone u ,afkcritical (λ − λ ) 3 These equations decouple in the harmonic approximation, i.e., 0 the TS t u s This has been computing velocity for systems described by anharmonic c 4 and second-order corrections, ass well asκ that of the higher ors˙ = s-3.2 . (21b) -4 term, the lead λu − λs c = − = k T . 4 B if f (x) = 0, 2but are still subject to the time-dependent u they s f (x) = 0, but they are still subject the demonstrated time-dependent 3 2 potentials using aif perturbative scheme. Furthermore, wetohave how can the critical dependent stochastic driving in -0.02 4k T ω λ 4 der terms that are omitted, varies among different realizations. ω λ (λ − λ ) B b u Downloaded 19 Jun 2012 to 158.125.33.105. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/about/rights_and_permissions s u s obtained from b be -6 stochastic driving force ξ (t). This time dependence can 1 α λ t velocity be used to calculate factor. this way, we have succeeded in obtaining s Nevertheless, the coordinate shif ¯ In -3.6 (56) -0.04 In the special case that the anharmonic potential contains "sk (t) S[λ = − "ukthe (0)etranmission − , f ; t], s k reads removed by -0.4 -0.2the coordinate 0 0.2 shift0.4 -8 λ −factor. λ analytical correction to the famous Kramers’

1

3

Since the critical velocity characterizes reactive trajectories, the transmission factor (7) can be expressed in terms of V ‡ (see Ref. 33 and Sec. V below). This picture of the invariant manifolds was introduced in Refs. 30 and 31 and applied to rate calculations in Refs. 32 and 33. The main purpose of the present work is to explore how this picture changes when anharmonicities of the barrier potential are taken into account. In this case the equations of motion (19) are coupled in a nonlinear time-dependent way, and they cannot be solved easily. However, as long as the coupling is sufficiently weak, one can expect to find a TS trajectory with its associated stable and unstable manifolds that are close to those in the harmonic approximation. Indeed, there are general theorems in the theory of stochastic dynamical systems42 that guarantee the persistence of these structures. As shown in Fig. 1(c), the invariant manifolds in an anharmonic system will beκtangent to their harmonic approximations at the TS trajectory, but they will not be straight lines anymore. Because the coupling term in Eqs. (19) is stochasti-

(23)

‡

V ‡ ≡ V0 = (λu − λs )u‡ (0).

!u

x

!s

(c)

κ

V‡

vx

!s

(x (t),v (t))

‡ ‡

B B

(x‡ (t),v‡ (t))

‡

‡

(x (t),v (t))

vx

‡ V Vcrit

FIG. 1. Phase space view of the time-dependent invariant manifolds of the Langevin equation. (a) Invariant manifolds are time-independent in the harmonic approximation andvyin relative coordinates. (b)vyIn space-fixed coordinates, the invariant manifolds are attached to the TS trajectory and move through phase space with it. (c) Anharmonic coupling deforms the manifolds. Both their position and their shape are stochastically time dependent. (d) Invariant manifolds can deviate strongly from the harmonic approximation if the anharmonicities are strong.

!s

(d)

x

!u

x

Downloaded 19 Jun 2012 to 158.125.33.105. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/about/rights_and_permissions

‡ Vcrit V

b

FIG. 1. Phase space view of the time-dependent invariant manifolds of the

(x‡ (t),v‡ (t))

(x (t),v (t))

‡ ‡

κ-κ0

vx

!s

(d)

x

!u

‡ Vcrit V

!u !x

non− reactive

b

4

b

3

B

-2.0

[email protected]

-2.4

!s

!s

κ2 3c k T =(37) − κ0 4 ωb4 V‡

vx

κ

0.56

1−µ 2 1 + µ -8 -6

(57)

-0.04

0.1

‡

0 c

‡

non− reactive

-0.1

(x (t),v (t))

A comparison -1 -0.5 0 0.5 1 1.5 2 2.5 3 in Fig. 5. Aga yλu /ωb that (c) in terms of the dimensionless parameter µ = κ = 0 (a) Critical velocity for one realization of the noise for Difference between numerically calculated critical velocity and If both c ω =1, ω =1.5, was γ=2, k used T=1. (b) in Transmission factor for Ref. 27. perturbative expansions for ω =1, ω =1.5, γ=2, c=0.2, k T=1. (a) ω =1, ω =0.5, γ=1, k T=1. the Harmonic approximation. (b) First-order in perturbationpotential, theory. comparison of Eq. (56) with numerical results is shown Numerical results (redAcrosses), harmonic Kramers) (c) Second-order perturbation theory. Contour spacing is 0.05 in approximation (gray horizontal line), perturbative in Fig. 4.straight They confirm that thein (b) perturbative result (a), and 0.005 and (c). Note that the color scaletransmission is also results to first –order (green line) and second-once more stretched by a factor 10 in (a). order (blue line). (57) and (59). is correct. The figure also shows the second-order correction 0.52 -0.2

-1.8

!u

(c)

x

!u !x

‡

4 2 -0.02 Downloaded 19 Jun 2012 to 158.125.33.105. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/about/rig 4 B

V2 = − c4 S[λu , X3 ; 0],

-2.2294974. *Financial support provided by MTM2012-39101-C02-01 and FP7/2007-2013 under Grant No. 3

f2 = − c4 X − 2c3 X "x1 .

x

Refs. 31. The equations of motion (19) decouple Since30theand critical velocity characterizes reactive trajectories, ‡ and the become time independent when f (x) = 0,inand they can transmission factor (7) can be expressed terms of V !uthen (see be easily by V writing Ref. 33solved and Sec. below). x in

s

(b)

(a)

reactive

which it canand be M. finally obtained [1] P.from Hänngi, P. Talkner, Borkovec, Rev. Mod. Phys. 62, 251 (1990). [2] P. Pechukas, in Dynamics ‡of Molecular Collisions, Part B, edited by W. H. Miller Vk pp. = 269-322; (λu − λD. (34) s )"u k (0). J. Chem. Phys. 68, (Plenum, New York, 1976) Chadler, 2959 (1978). [3] W. H. Miller, Faraday Discuss. Chem. Soc. 110, 1 (1998). recursion relations canPhys. be successively evaluated as [4] T.The Bartsch, R. Hernandez, and(33) T. Uzer, Rev. Lett. 95 058301 (2005). = 1, 2, . and . . up to any desired order. [5] T.written Bartsch,for R. kHernandez, T. Uzer, J. Phys. Chem. 123, 204102 (2005). [6] T. Bartsch, T. Uzer, J. M. and R. Hernandez, Chem. B 112, For example, forMoix, the anharmonic forceJ. Phys. corresponding to206 (2008). [7] F.the Revuelta, Bartsch, R. M. Benito, F. Borondo, J. Chem. genericT. one-dimensional potential (10) with Phys. only (Communication) cubic 136, 091102 (2012). and quartic terms, expansion (32) gives [8] T. Bartsch, F. Revuelta, R. M. Benito, F. Borondo, J. Chem. Phys. 136, 224510 (2012). [9] H. A. Kramers, Physica = − c X7,2 ,284 (1940). f (Ultrecht)

!v

REFERENCES

(33)

4

driving from the leading-order t c only a quartic term, the perturbation expansion results as an 0.68 -1 -0.5 0 0.5 1 1.5 2 2.5 3 0.04 0 ‡ ‡ (c) it into the anharmonic perturbati 2 , (s = s − s , (16) (u = u − u y (b) odd orders terms null. For expansion in powers of ε , with theThis (b) result agrees with the perturbative correction given in c3 0.64 0.02 -2 c κ 1 the first two non-zero corrections, aRefs. similar calculation shows 2 26–28. It can be rewritten as =− 0 0.60 -4 that κ0 6 ! "2 c

Vcrit

s

4

3

3

vx

Since λu > 0 and λs < 0, the coordinate !u grows exponentially in time, whereas !s shrinks. Therefore, !u and !s correspond to unstable and stable directions in phase space, respectively. In particular, the lines !u = 0 and !s = 0 are invariant under the dynamics. A trajectory that starts on the line !u = 0 will asymptotically approach the origin as t → ∞; this line is called the stable manifold of the origin. A trajectory on the line !s = 0 will move away from the origin as t → ∞, but it will approach the origin as t → −∞; this line is called the unstable manifold of the origin. The stable and unstable manifolds of the origin, together with several typical trajectories in relative coordinates, are shown in Fig. 1(a). The invariant manifolds separate trajectories with different qualitative behavior. Trajectories above the stable manifold, i.e., with larger relative velocity, move to the product side of the barrier for asymptotically long times, whereas trajectories below the stable manifold move to the reactant side. Similarly, trajectories above the unstable manifold come from the reactant side in the distant past, whereas

(22)

!s(t) = !s(0) eλs t .

!u(t) = !u(0) eλu t ,

The description of the geometric phase space structure in the vicinity of the saddle point is most easily done if one starts from the harmonic limit. A full discussion can be found in Refs. 30 and 31. The equations of motion (19) decouple and become time independent when f (x) = 0, and they can then be easily solved by writing

reactive ‡

‡

4

4

= (λu

‡

u

4

V ≡

x

‡

u

!u − λ!x)u‡ (0).

‡

‡ V0

‡

"xk (t) = "uk (t) + "sk (t), Bartsch et al.

(23)

‡

u

Since the critical velocity characterizes reactive trajectories, ‡ ‡ !s (x (t),v (t)) the transmission factor (7) can be expressed in terms of V ‡ This picture !u(t) of the invariant was introduced in = !u(0)manifolds e λu t , (see Ref. 33 and Sec. V below). Refs. 30 and 31 and applied to rate calculations in Refs. 32 = !s(0) eλs t . work is to explore (22) and 33. The main!s(t) purpose of the present This picture of the invariant manifolds was introduced in!s when anharmonicities of the barrier Refs. 30 and 31 and applied to rate calculations in Refs. 32 Sincehow λu this > 0picture and λschanges < 0, the coordinate !u grows exponenarewhereas taken into!s account. In this case the equations of tiallypotential in time, shrinks. Therefore, !u and !s and 33. The main purpose of the present work is to explore vx vx barrier motion (19) are coupled a nonlinear time-dependent way, to unstable and instable directions in phase space, how this picture changes when anharmonicities of the (c) (d) correspond respectively. In particular, lines !u =as0long andas !s = 0 and they cannot be solved the easily. However, the coupotential are taken into account. In this case the equations of are invariant under theweak, dynamics. trajectory starts on pling is sufficiently one canAexpect to findthat a TS trajecmotion (19) are coupled in a nonlinear time-dependent way, the line will asymptotically approach the origin as tory !u with= its 0associated stable and unstable manifolds that are and they cannot be solved easily. However, as long as the cout → close ∞; this line isincalled the stable manifold of Indeed, the origin. to those the harmonic approximation. thereA x x pling is sufficiently weak, one can expect to find a TS trajectrajectory on thetheorems line !s in =0 move fromdynamical the origin are general thewill theory of away stochastic ‡ tory with its associated stable and unstable manifolds that are V as t → ∞, 42 butthat it will approach the originofasthese t → structures. −∞; this ‡ ‡ !u guarantee the persistence systems (x (t),v (t)) ‡ ‡ (t),v harmonic (t)) close to those in(xthe approximation. Indeed, there line is called the unstable manifold of the origin. in an anharAs shown in Fig. 1(c), the invariant manifolds κ The stable andwill unstable manifolds ofharmonic the origin, together are general theorems in the theory of stochastic dynamical monic system be tangent to their approximawith several typical trajectories in relative coordinates, are !s of these structures. systems42 that guarantee the persistence tions at the TS trajectory, but they will not be straight lines shown in Fig. 1(a). The invariant manifolds separate trajecAs shown in Fig. 1(c), the invariant manifolds in an anharanymore. Because the coupling term in Eqs. (19) is stochastiFIG. 1. Phase space view of the time-dependent invariant manifolds of thetories with different qualitative behavior. Trajectories above monic system Langevin will beequation. tangent to their harmonic approximacally time dependent, the shapes of the invariant manifolds as (a) Invariant manifolds are time-independent in the har-the stable manifold, i.e., with larger relative velocity, move to tions at the TSmonic trajectory, but they will not be straight lines well as their positions in phase space depend on time and on approximation and in relative coordinates. (b) In space-fixed coordithe product side of the barrier for asymptotically long times, nates, the the invariant manifolds to the is TS stochastitrajectory and move anymore. Because coupling termareinattached Eqs. (19) the realization of the noise.the stable manifold move to the whereas trajectories below through phase space with it. (c) Anharmonic coupling deforms the manifolds. The intersection of the stable manifold with the axis cally time dependent, the shapes of the invariant manifolds as Both their position and their shape are stochastically time dependent. (d) In-reactant side. Similarly, trajectories above the ‡unstable manisuch that tra= 0 from will give to a critical V past, in phase space depend onharmonic time and on vywell as their positions variant manifolds can deviate strongly from the approximation iffold xcome the rise reactant side in velocity the distant whereas jectories with initial velocities larger than V ‡ will be reactive, the realization the ofanharmonicities the noise. are strong. The intersection of the stable manifold with the axis vx !v ‡ (a) x = 0 will give rise to a critical velocity V such that tra(b) ‡ Downloaded 19 Jun 2012 to http://jcp.aip.org/about/rights_and_permissions will be Redistribution reactive,subject to AIP license or copyright; seereactive jectories with initial velocities larger than V158.125.33.105.

trajectories below the unstable manifold come from the product side. For a reaction rate calculation we need to ascertain whether a trajectory will turn into reactants or products in the future. In our approach this sentence is rephrased into the condition: We need to decide whether a trajectory lies above or below the stable manifold. In other words, the stable manifold encodes the information about the reaction dynamics that is most relevant to us. We will therefore focus on the stable manifold in what follows, largely ignoring the unstable manifold. We can return to space fixed coordinates by undoing the time-dependent shift (16). After the shift, the stable and unstable manifolds are not attached to the origin of the coordinate system any more, but instead to the TS trajectory as a moving origin, as shown in Fig. 1(b). Since the TS trajectory is time dependent, the manifolds will move through phase space with it. Nevertheless, they still separate trajectories with different asymptotic behaviors. Given a trajectory with a given initial condition at a certain time, it can be classified as reactive or non-reactive by knowing the instantaneous position of the stable manifold at that time. Through the TS trajectory, that instantaneous position will depend on the realization of the noise. It is clear from Fig. 1(b) that at any time and for any realization of the noise the stable manifold intersects the axis x = 0 at a point with a velocity V ‡ . Trajectories with initial positions x = 0 and initial velocities vx > V ‡ are reactive, while trajectories with initial velocities vx < V ‡ are not. The critical velocity V ‡ depends on time and on the realization of the noise. For the harmonic approximation, it was shown in Ref. 33, and it will be rederived below, that

J. Chem. Phys. 136, 224510 (2012)

Vcrit

trajectories belowThe thedescription unstableof manifold come the geometric phasefrom space the structure trajectories below the unstable manifold come from the product side. in the vicinity of the saddle point is most easily done if one product side. For a reaction ratethecalculation weA full need to ascertain starts from harmonic limit. discussion can be found For a reaction rate calculation we need to ascertain whether a trajectory turn reactantsofor products in in Refs. will 30 and 31. into The equations motion (19) decouple whether a trajectory will turn into reactants or products in the future. In our thisindependent sentence is rephrased and approach become time when f (x) = 0,into and the they can the future. In our approach this sentence is rephrased into the solved by writing condition: We need to decide whether a trajectory lies above condition: We then needbetoeasily decide whether a trajectory lies above or below the stable manifold. In other words, the stable manior below the stable manifold. In!u(t) other=words, λu t stable mani, !u(0) ethe fold encodes the information about the reaction dynamics that fold encodes the information about the reaction dynamics that λs t !s(t) = !s(0) e . on the sta- (22) is most relevant to us. We will therefore focus on the stais most relevant to us. We will therefore focus ble manifold in what follows, largely ignoring the unstable ble manifold in what follows, unstable Since λu > 0 and λs largely < 0, the ignoring coordinate the !u grows exponenmanifold. manifold. tially in time, whereas !s shrinks. Therefore, !u and !s We can return to space fixed coordinates by undoing the We can return to space fixed coordinates by undoing thespace, correspond to unstable and stable directions in phase time-dependent shift (16). After the shift, the stable and untime-dependentrespectively. shift (16).InAfter the shift, the !u stable particular, the lines = 0and andun!s = 0 stable manifolds are not attached to the origin of the coorare are invariant the dynamics. A trajectory starts on stable manifolds not under attached to the origin of thethat coordinate system any more, but instead to the TS trajectory as the line !u = 0 will asymptotically approach the origin as dinate system any more, but instead to the TS trajectory as a moving origin, as shown in Fig. 1(b). Since the TS traject → ∞; this line is called the stable manifold of the origin. A a moving origin, as shown in Fig. 1(b). Since the TS trajectory is time dependent, the manifolds will move through phase trajectory on the line !s = 0 will move away from the origin tory is time dependent, the manifolds will move through phase space with it. Nevertheless, they still separate trajectories as t → ∞, but it will approach the origin as t → −∞; this space with it. Nevertheless, they still separate trajectories with different asymptotic behaviors. Given a trajectory with a line is called the unstable manifold of the origin. with different asymptotic behaviors. Given a trajectory with a given initial condition at a certain time, it can be classified as The stable and unstable manifolds of the origin, together given initial condition at a certain time, it can be classified as reactive or non-reactive by knowing the instantaneous posiwith several typical trajectories in relative coordinates, are reactive or non-reactive by knowing the instantaneous position of the stable manifold at that time. Through the TS trajecshown in Fig. 1(a). The invariant manifolds separate trajectory, that instantaneous position will depend on the realization tion of the stable manifold at that time. Through the TS trajectories with different qualitative behavior. Trajectories above of the noise. tory, that instantaneous position will depend on the realization the stable manifold, i.e., with larger relative velocity, move to It is clear from Fig. 1(b) that at any time and for any of the noise. the product side of the barrier for asymptotically long times, realization of the noise the stable manifold intersects the axis It is clearwhereas from Fig. 1(b) that anystable timemanifold and formove anyto the trajectories belowat the ‡ . Trajectories with initial x = 0 at a point with a velocity V realization of the noise theSimilarly, stable manifold axismanireactant side. trajectoriesintersects above the the unstable ‡ > V are reactive, positions x = 0 and initial velocities v ‡ x withpast, initial x = 0 at a point velocity V . Trajectories foldwith comea from the reactant side in the distant whereas while trajectories with initial velocities vx < V ‡ are not. The ‡ positions x = 0 and initial velocities vx > V are reactive, ‡ depends on time and on the realization of critical velocity 224510-5 BartschVet al. ‡ while trajectories with initial velocities vx < V are not. The vx the noise. For the harmonic approximation, it was shown in !v ‡ critical velocity V depends on time(a)and on the realization of (b) Ref. and it willof be the rederived below,phase that space structure The 33, description geometric the noise. For the harmonic reactive approximation, it was shown in in the vicinity of the ‡saddle‡ point is most easily done if one (23) V ≡ V0 = (λu − λs )u‡ (0). Ref. 33, and it will be rederived below, that non− starts from the harmonic limit. A full discussion can be found

J. Chem. Phys. 136, 224510 (2012) 224510-5

Bartsch et al.

J. Chem. Phys. 136, 224510 (2012)

x

224510-5

1 2 2 β 4 4 x y + (x + y ) 2 4

Reaction rate calculation for dissipative systems using invariant manifolds

x

y

x

y

B

B

0.2

x

y

B