Transition State Theory for dissipative systems without a dividing surface
Transition State Theory is a central cornerstone in reaction dynamics. Its key step is the identification of a dividing surface that is crossed only once by all reactive trajectories. This assumption is often badly violated, especially when the reactive system is coupled to an environment. The calculations made in this way then overestimate the reaction rate and the results depend critically on the choice of the dividing surface. In this Letter, we study the phase space of a stochastically driven system close to an energetic barrier in order to identify the geometric structure unambiguously determining the reactive trajectories, which is then incorporated in a simple rate formula for reactions in condensed phase that is both independent of the dividing surface and exact.
pacs:82.20.Db, 05.40.Ca, 05.45.-a, 34.10.+x
Transition State Theory (TST) provides the conceptual framework for a large part of reaction rate theory. Originally developed to describe reactivity of small molecules Truhlar, Hase, and Hynes (1983); Truhlar, Garrett, and Klippenstein (1996); Miller (1998), it was later extended to study a wide variety of processes in very different fields that only have in common the existence of a transition from well-defined “reactant” to “product” states Toller et al. (1985); Eckhardt (1995); Jaffé, Farrelly, and Uzer (2000); Koon et al. (2000); Jaffé et al. (2002); Uzer et al. (2002); Komatsuzaki and Berry (2002); Waalkens, Burbanks, and Wiggins (2004a, b). Its great success comes from its simplicity, since it gives a straightforward answer to the two central problems in chemical dynamics: the identification of the reaction mechanism and a simple approximation to the reaction rate.
The rate limiting step in many reactions is the crossing of an energetic barrier, the top of which forms a bottleneck that the system must cross as the reaction takes place. If a dividing surface (DS) is placed close to this bottleneck, reaction rates can be readily computed from the steady-state flux through it. The TST approximation is obtained under the assumption that the reactive classical trajectories cross the DS only once and never return. Very often, for example when the system is strongly coupled to a noisy environment such as a liquid solvent, this no-recrossing assumption fails, and any conventional DS is crossed many times by a typical trajectory. As a result, TST calculations significantly overestimate reaction rates, and considerable effort has been devoted to the construction of a DS that minimizes recrossing Truhlar, Garrett, and Klippenstein (1996). The problem mainly derives from the fact that this special surface plays a double role: It defines (or separates) reactant and product regions, and it also serves to identify the reactive trajectories. To the latter purpose the DS is ill suited. Despite this fundamental drawback, TST remains attractive since reactive trajectories are simply identified as those crossing the DS. This criterion only takes into account the instantaneous velocity at the DS, without requiring any time consuming trajectory simulation.
Important advances in TST have recently been achieved within the approach of modern nonlinear dynamics: A strictly recrossing free DS can be constructed in the phase space of reactive systems with arbitrarily many degrees of freedom Uzer et al. (2002); Waalkens, Burbanks, and Wiggins (2004a, b). These results were generalized to systems interacting with an environment Bartsch, Hernandez, and Uzer (2005); Bartsch, Uzer, and Hernandez (2005); Bartsch et al. (2006, 2008); Hernandez, Uzer, and Bartsch (2010); Kawai and Komatsuzaki (2010a, b). It has been shown that the desired exact TST can be constructed, in the harmonic limit, by using a moving DS only crossed once by the reactive trajectories Bartsch, Hernandez, and Uzer (2005); Bartsch, Uzer, and Hernandez (2005). Accurate results are still obtained for moderate anharmonicities Bartsch et al. (2006) and can be improved by normal form calculations Kawai and Komatsuzaki (2010a, b).
In this Letter, we make a step forward by presenting an explicit description of geometric phase space structures in an anharmonic noisy system and an analytical scheme that relies on these structures for the computation of exact TST reaction rates for arbitrary multidimensional potentials coupled to a noisy environment. The key point is the demonstration that reactive trajectories can be rigorously identified solely from their initial conditions, thus avoiding the choice of a (arbitrary) DS. This is done in terms of the stable manifold associated to a “noisy” Transition State (TS) trajectory jiggling in phase space. This geometric structure encodes the relevant information about the noise in the most economical manner and it can easily be incorporated into a rate calculation. Our method retains the fundamental simplicity of TST, also providing the conceptual tools to develop new computational algorithms. An application to the simple case of the one-dimensional quartic potential is presented as an illustration.
The Langevin equation (LE) has been widely used to model the interaction of a reactive system with a surrounding heat bath Hänggi, Talkner, and Borkovec (1990); Pechukas (1976); Chandler (1978). Being a classical model, this description neglects quantum effects such as barrier tunnelling, which can be important in the case of light particles Bothma, Gilmore, and McKenzie (2010), and the interaction with excited surfaces through conical intersections Polli and et al. (2010).
The dynamics of a unit mass particle defined by coordinate moving in a one-dimensional potential is given by
Here, is the fluctuating force exerted by the bath, which is connected to the damping strength, , by the fluctuation–dissipation theorem
The model potential that we have chosen to study is
although our derivation equally applies to any other one-dimensional case. This generality will be emphasized by using , equal to in our case, to denote any anharmonic force. The extension to higher dimension is straightforward and will be presented elsewhere.
For every fixed realization of the noise, the LE gives rise to a specific trajectory called TS trajectory Bartsch, Hernandez, and Uzer (2005); Bartsch, Uzer, and Hernandez (2005). This orbit remains in the vicinity of the energetic barrier for all times, without ever descending into any of the potential wells. Other trajectories will be referred to the TS trajectory, which will be taken as a moving coordinate origin. Since the LE is a second order differential equation, its phase space is two-dimensional, with coordinates and . If we now introduce the new coordinates
relative values are defined by the time-dependent shift
where the coordinates of the TS-trajectory are
The corresponding equations of motion are
with =+ [and =+]. Now, the geometric phase space structure in the vicinity of the barrier can be easily discussed. First, in the harmonic limit, i.e. , the equations of motion (8) are trivially solved, giving
Since and , increases exponentially in time, whereas shrinks accordingly. More importantly, the lines and are invariant under the dynamics of the system; being the unstable and stable manifolds of the origin, respectively. As shown in Fig. 1(a), these invariant manifolds separate trajectories with different qualitative behavior: Those above the stable manifold (larger relative velocity) move to the product side in the distant future, while those below it move, on the other hand, towards the reactant side.
In space fixed coordinates the invariant manifolds appear attached to the TS trajectory, as shown by the dashed lines in Fig. 1(b); their instantaneous position depends on the realization of the noise. Accordingly, the manifolds move through phase space but they still separate trajectories with different asymptotic behaviors. The stable manifold intersects the axis in a point with velocity . Trajectories with initial position and initial velocity larger than this critical velocity, , are reactive, while trajectories with initial velocities are not. The (random) value therefore encodes the relevant information about the realization of the noise concisely. In other words, once the instantaneous position of the stable manifold is known, any trajectory can unambiguously be classified from the values of its initial condition as reactive or non-reactive. Finally, the presence of anharmonicities () will distort the invariant manifolds, as indicated by the red lines with arrows in Fig. 1(b). The main step of the theory to be developed here is the calculation of this deformation.
This critical velocity can be calculated from the condition that the trajectory with and is contained in the stable manifold of the TS trajectory. To find this trajectory, which will be called critical trajectory, we formally solve the equations of motion (8) by
where and are two arbitrary constants, and the integral operator
has been introduced as a convenient shorthand notation. The subscript indicating the integration variable in Eq. (11) will be left out whenever this does not cause any ambiguity. The unknown constants and can be determined by noticing that the critical trajectory should approach the TS trajectory for large times. In particular, it should remain bounded for . This can only be satisfied if . [This condition also ensures that the functional in Eq. (10a) is well defined.] With this choice, Eq. (10a) determines the initial condition , and can then be found from the condition that . Finally, the coordinate transformation (4) yields
In general, Eq. (10) represents only a formal solution, since its right-hand side depends on the unknown functions and . However, since it is known that the critical trajectory remains in the neighborhood of the barrier at all times, the anharmonic force will be small, and then its influence can be evaluated through perturbation theory, thus obtaining an expansion in powers of the anharmonic coupling parameter .
In the harmonic approximation the critical trajectory is given by
and for this case Eq. (12) yields
which is the harmonic approximation to the coordinate of the critical trajectory. Equation (10a) then gives
and therefore the leading-order velocity correction is
is obtained. In Fig. 2 we present a comparison between numerical (exact) results for the critical velocity, , for one realization of the noise and the approximate values computed using (15) and (16), as a function of the anharmonic parameter, .
In order to calculate the corresponding reaction rate, , we choose the simplest DS, defined by , and use the basic flux-over-population rate formula (see e.g. Hänggi, Talkner, and Borkovec (1990)), which states that is proportional to the reactive flux
across the DS. This flux is to be averaged over different realizations of the noise, , and also over a Boltzmann ensemble of initial velocities, , for trajectories starting at the DS. Notice that only reactive trajectories should be included in the average.
The non-recrossing assumption of conventional TST can be restated here by saying that the reactive trajectories are those crossing the DS with velocity . We call the rate constant obtained with this approximation . Any effects beyond TST are customarily summarized Hänggi, Talkner, and Borkovec (1990) into a transmission coefficient .
In terms of our stochastic invariant manifolds, reactive trajectories are characterized by , as discussed before. Using this criterion, the Boltzmann average over velocities in Eq. (17) can be evaluated, as it was in Ref. Bartsch et al., 2008 for the harmonic case. This gives the exact expression
where only the average over the noise remains to be done. By substituting the perturbative expansion for the critical velocity into this expression and expanding the exponential, a perturbative series of rate corrections, , is obtained, where
with the abbreviated notation
Expressions similar to (19) can be obtained for the higher-order corrections. The leading order was evaluated in Ref. Bartsch et al., 2008. It yields the well-known Kramers result for the transmission factor. With the result (15), the first correction term is given by
To perform the remaining average and subsequently evaluate the functional, note that and , and consequently , are Gaussian random processes, whose correlation functions were given in Bartsch, Uzer, and Hernandez, 2005. Details of this calculation are irrelevant for the purpose of this Letter and will be presented elsewhere. The final result is given by
in terms of the dimensionless parameter . This expression agrees with the corrections given in Refs. Pollak and Talkner, 1993; Talkner and Pollak, 1993; Talkner, 1994. A comparison between perturbation theory and the results of a numerical simulation is shown in Fig. 3, along with the second-order perturbative correction, which can be obtained in a similar manner. The perturbative results describe the rate correctly as long as the coupling is not too strong. For large negative values of the second-order correction loses its accuracy. By contrast, it is accurate for all positive shown in the figure. If is increased further, the wells of the model potential (3) become too shallow for a rate theory to be meaningful.
In summary, we have demonstrated that the use of stochastic invariant structures allows a rigorous TST reaction rate calculation in general anharmonic and noisy systems. In our approach, the DS is only used to define reactant and product regions. To identify reactive trajectories we employ the stable manifold that is determined by the dynamics of the system itself. The resulting rate formula (18) is not only remarkably compact, but also exact. In particular, the arbitrariness that is usually introduced by the choice of a DS is absent. Although our presentation is based on an analytical perturbation expansion, the method can easily be incorporated into a numerical scheme to achieve an efficient rate calculation in complex systems.
Support from MICINN–Spain under Contracts nrs. MTM2009–14621 and i–MATH CSD2006–32 is gratefully acknowledged. FR thanks UPM for a doctoral fellowship and the hospitality of the members of the School of Mathematics at Loughborough University, where part of this work was done.
- D. G. Truhlar, W. L. Hase, and J. T. Hynes, J. Phys. Chem. 87, 2664 (1983).
- D. G. Truhlar, B. C. Garrett, and S. J. Klippenstein, J. Phys. Chem. 100, 12771 (1996).
- W. H. Miller, Faraday Discuss. Chem. Soc. 110, 1 (1998).
- M. Toller, G. Jacucci, G. DeLorenzi, and C. P. Flynn, Phys. Rev. B 32, 2082 (1985).
- B. Eckhardt, J. Phys. A 28, 3469 (1995).
- C. Jaffé, D. Farrelly, and T. Uzer, Phys. Rev. Lett. 84, 610 (2000).
- W. S. Koon, M. W. Lo, J. E. Marsden, and S. D. Ross, Chaos 10, 427 (2000).
- C. Jaffé, S. D. Ross, M. W. Lo, J. Marsden, D. Farrelly, and T. Uzer, Phys. Rev. Lett. 89, 011101 (2002).
- T. Uzer, C. Jaffé, J. Palacián, P. Yanguas, and S. Wiggins, Nonlinearity 15, 957 (2002).
- T. Komatsuzaki and R. S. Berry, Adv. Chem. Phys. 123, 79 (2002).
- H. Waalkens, A. Burbanks, and S. Wiggins, J. Phys. A 37, L257 (2004a).
- H. Waalkens, A. Burbanks, and S. Wiggins, J. Chem. Phys. 121, 6207 (2004b).
- T. Bartsch, R. Hernandez, and T. Uzer, Phys. Rev. Lett. 95, 058301 (2005).
- T. Bartsch, T. Uzer, and R. Hernandez, J. Chem. Phys. 123, 204102 (2005).
- T. Bartsch, T. Uzer, J. M. Moix, and R. Hernandez, J. Chem. Phys. 124, 244310 (2006).
- T. Bartsch, T. Uzer, J. M. Moix, and R. Hernandez, J. Phys. Chem. B 112, 206 (2008).
- R. Hernandez, T. Uzer, and T. Bartsch, Chem. Phys. 370, 270 (2010).
- S. Kawai and T. Komatsuzaki, PCCP 12, 7626 (2010a).
- S. Kawai and T. Komatsuzaki, PCCP 12, 7636 (2010b).
- P. Hänggi, P. Talkner, and M. Borkovec, Rev. Mod. Phys. 62, 251 (1990).
- P. Pechukas, in Dynamics of Molecular Collisions, Part B, edited by W. H. Miller (Plenum, New York, 1976).
- D. Chandler, J. Chem. Phys. 68, 2959 (1978).
- J. P. Bothma, J. B. Gilmore, and R. H. McKenzie, New J. Phys. 12, 055002 (2010).
- D. Polli et al., Nature 467, 440 (2010).
- E. Pollak and P. Talkner, Phys. Rev. E 47, 922 (1993).
- P. Talkner and E. Pollak, Phys. Rev. E 47, R21 (1993).
- P. Talkner, Chem. Phys. 180, 199 (1994).