Analytic dynamical corrections to transition state theory

We present a theory of reaction kinetics that accurately accounts for recrossing effect at the transition state. Using perturbation theory we solve the equations of motion in the neighborhood of the transition state and derive a simple condition that determines whether a trajectory is truly reactive or not. We apply our theory to adatom diffusion on an aluminum surface and show that our analytical results compare well with those of a direct dynamical simulation. Our theory gives us a better understanding of the limitations of transition state theory and a simple analytic formula for quantifying the recrossing effects.

Transition state theory (TST) plays a crucial role for understanding kinetic phenomena in chemistry and physics, [1][2][3] including reaction kinetics, diffusion, and transport. TST is also essential for simulating long time-scale dynamics, [4][5][6] where it can be used to avoid expensive molecular dynamics (MD) calculations for slow processes. Notably, the harmonic approximation of TST has proven to be extremely useful because of its simplicity. Within harmonic TST (HTST) [7,8] one only needs to determine the barrier height for a reaction and the vibrational frequencies for the initial state and the transition state (TS). Finding new methods for improving TST or HTST can significantly improve efforts to make realistic predictions of kinetics for complex and large systems. The current paper makes a step toward refining TST by considering contributions from stochastic interactions at the TS due to anharmonic terms in the Hamiltonian.
The fundamental assumption of TST is that when a trajectory crosses the TS, there are no subsequent recrossings. Therefore, TST variationally overestimates the true rate. To see this in detail, we can look at a general expression for the reaction rate where θ is the Heaviside step function, x is the reaction coordinate, and the TS is taken to be at x = 0 (see figure 1). For rare events, the characteristic time τ in this equation is much smaller than the reaction time, ∼ k 1 . For simplicity of notation, when we write the position and velocity without a time argument, we mean that they are the initial position and velocity at time zero. The angular brackets indicate an ensemble average at time zero. The denominator in equation (1) is the partition function of the reactant state.
The reason that equation (1) gives a more accurate rate than TST is that it does not over-count recrossing events when calculating the flux through the TS. Consider a trajectory that starts at x = 0 with positive velocity which is scattered back across the TS (see figure 1). TST includes this recrossing as contributing to the TST rate, whereas equation (1) does not include this unsuccessful reaction in the true rate. Consider a second trajectory that also starts at x = 0 with positive velocity, but is coming from the product side rather than the reactant side. This is obviously not a true reactive crossing as it comes from the product state and goes back to the product state. Again, TST treats this as a reactive crossing. This trajectory also contributes to equation (1) with positive sign. However, a similar trajectory that starts on the product side, crosses the TS, and moves back to the product, with a negative velocity at time zero, contributes negatively to the rate. These positive and negative contributions cancel out so that there is no net contribution to the true rate from recrossing events.
Many authors have worked to generalize the concept of a dividing surface to address the recrossing problem. One of the most direct ways to define an optimal dividing surface is based on variational transition state theory (VTST) [9,10]. VTST utilizes the fact that TST gives an upper limit of the rate, so that within VTST, the dividing surface is defined as a surface across which the reaction flux is minimal. Simulations on realistic systems, however, show that VTST cannot fully eliminate recrossing effects [26,27]. One of the earliest studies that provides a microscopic understanding of dynamics near the saddle point was done by Pollak and coworkers [11,12]. Reference [13] provides an excellent review of early work in this area. Later Jacucci and coworkers derived an analytic expression for the transmission coefficient for non-planar dividing surfaces [28,29,41]. More recently, Uzer and Hernandez defined a rigorous mathematical approach to extend the definition of a dividing surface beyond configuration space to phase space [13][14][15][16][17][18]. Defined in phase space, these authors show that a dividing surface can be exact for some model potentials. An alternative approach to calculating a transmission coefficient is Grote-Hynes theory [30] and extensions [31] based upon the generalized Langevin equation (GLE) [19,32]. While Grote-Hynes theory has been used extensively, [33,34] it requires the computation of a memory kernel. Other important developments related to the TST dividing surface and transmission coefficient can be found in [35][36][37][38].
The main idea behind our approach is similar to that of the Uzer-Hernandez method. Our calculation of a transition state trajectory has also been obtained by other authors [28]. Our technical implementation, however, is different from previous work and is more suitable to our problem. Specifically, our objective is to obtain an analytic expression for the transmission coefficient in terms of force constants. Such an expression eliminates the need for calculating molecular dynamics trajectories or any other expensive simulation.
The paper is organized as follows. First we introduce expressions for the transmission coefficient κ. Then, using perturbation theory we study trajectories near the TS and use a mathematical analysis of these trajectories to obtain an analytic expression for κ. Finally, we validate our theory on a two-dimensional potential and a higher dimensional system involving atomic diffusion on a metal surface.

Transmission coefficient
The transmission coefficient κ was introduced to describe deviations of the TST rate from the true rate of a reaction q t d q d Here, we have multiplied and divided equation (1) by the flux through the TS at time zero. This gives us the definition of κ Equation (3) is the primary equation that we aim to approximate in this work. While a numerical approach for determining κ has been already implemented, [8,20] our goal is to find an analytic solution for κ and also to obtain a better understanding of the stochastic processes that lead to recrossing events in TST.

Anharmonic motion near the TS
Since κ is a property of the TS, we will focus on trajectories in the vicinity of the TS. The total Hamiltonian of the system can be split into harmonic and anharmonic terms. The harmonic Hamiltonian, in mass weighted normal coordinates, is where Q 0 corresponds to the negative mode at the saddle point, and the remaining modes from i = 1 to N are positive. Q and P are related to the real positions and momenta, q and p, by where R is the matrix that diagonalizes the Hessian matrix Next we add anharmonicity, considering only terms of the type Q Q Q i j 0 . It can be shown that, for the perturbation theory we consider in this paper, terms involving higher orders of Q 0 or products of three positive modes (Q Q Q i j k ) do not have a significant effect on the trajectories that we are interested in. Ignoring the higher order terms, the interaction Hamiltonian is The coupling constants L ij can be expressed in terms of the derivatives of the Hessian matrix with respect to the negative mode along the reaction coordinate To find Λ, we displace the system along the negative mode to determine changes of the Hessian matrix within a finite difference approximation.
Under this anharmonic Hamiltonian, the equations of motion for Q 0 and Q i are We are interested in solving the first equation for ( ) Q t 0 . If we solve the second equation first however, and determine Q i (t), the second term in equation (10) can be written as an external force acting on Q 0 . This can be seen by explicitly writing equation (10) (see the appendix) as In other words, we are using perturbation theory and keeping only first order L ij terms in equation (12). For convenience of notation, we have introduced generalized position-momentum variables We will focus only on the first exponentially diverging term and avoid writing the other terms explicitly. In equation (14), G is aŃ where A, B, and C are N×N matrices Equation (14) tells us which trajectories successfully cross the TS and which will recross. If the factor -P G X X T 0 of the exponentially diverging term is positive, Q 0 increases to infinity, and the system successfully crosses the TS. The condition for the successful trajectory is If the factor is negative, Q 0 goes to negative infinity, and the crossing is unsuccessful. Over time w e t 0 goes to zero and the last term describes oscillatory motion around the TS. Additionally, if - , then once the system is near the TS, it will remain in that neighborhood. This critical trajectory that oscillates at the TS was introduced in [16] as a TS trajectory. Whether the trajectory is a TS trajectory or not depends on the initial positions and momenta. If the initial conditions change very slightly, the trajectory will no longer be a TS trajectory, and it will either go to the reactant or product states. Note that for the TS trajectory, perturbation theory is accurate. This is because Q 0 and the Q i ʼs approach zero near to the TS and the interaction Hamiltonian (equation (12)) is small at all times. As the system moves away from the TS, the perturbation expansion becomes inaccurate, but we are not interested what happens far from the TS because our assumption is that recrossing events can only occur in the vicinity of the TS. If the system moves away from the TS, falls into another basin, and after many random collisions recrosses the TS, its motion will be uncorrelated with the initial conditions and can be considered an independent trajectory.

Ensemble average
To obtain the transmission coefficient we require a canonical average over all trajectories that cross the TS. For this, we choose Q 0 as the reaction coordinate. Since we know which trajectories go to positive and negative infinity, we can rewrite equation (3) as We assume that the Hamiltonian is harmonic in the TS hyper-plane. This makes integration of the denominator of equation (20) straightforward since it is a Gaussian integral. To simplify the numerator, we can use an integral representation of the step-function Performing the integration with P 0 gives, where we used the approximation that 1, as there is no need to keep terms that contain infinitesimally small ò. In the exponent we have quadratic functions plus small correction terms. This type of integrals can evaluated using standard many body physics methods [21], such as the Feynman diagrams shown in figure 2. Although our system is classical, the formalism is very similar to that of quantum case. Our discussion closely follows that of the 'stationary phase approximation' in [21].
In the many body picture, is convenient to introduce new variables ... 0 indicates a Gaussian average. The auxiliary variable η in this equation behaves like another degree of freedom that interacts with all other degrees of freedom. Since b G is a small number, we can use it as a parameter in the perturbative expansion. To perform the integration we use Wick's theorem (known also as Isserlis' theorem in probability theory) and the linked cluster theorem (LCT). If we expand the exponent in equation (26) in a power series, the first order term is zero. We can then apply Wick's theorem up to the second order term  figure 2(a)). The former is equal to d ij , while the later is unity. Each vertex is connected to one dashed line and two solid lines. Each vertex gives a b -G i ij factor. We also need to multiply each diagram by ! S n , where S is the symmetry factor; the number of ways we can draw each diagram, and n is number of vertices. Feynman diagrams corresponding to equation (28) are shown in figure 2(b). We adopt the same terminology as in quantum many body physics and refer to the first one as the Hartree diagram and the second as the Fock diagram. If we think of the dashed line as an interaction between two points, than the analogy with the Hartree and Fock energies are clear. The Hartree interaction has symmetry factor one and the Fock interaction has two, which also can be seen from equation (28), where the two Fock terms are identical. The LCT avoids disconnected diagrams that appear in the Taylor expansion of the exponent in equation (26). According to the LCT, the logarithm of κ is the sum of all connected diagrams. Hence, within the Hartree-Fock (HF) approximation κ is k = Since the HF is a low order expansion, it can fail at high temperature. The simplest way to add higher order corrections is by using a Dyson equation for the dashed line. We compute the screened hh á ñ s correlation function using the simplest 'polarization function' shown in figure 2 The physical interpretation of screening is that interaction between two vibrational modes excites other vibrations and this reduces interaction strength. Since there is only one dashed line in both the Hartree and Fock diagrams, we obtain The subscript RPA stands for the random phase approximation [22,23]. Screening terms could also be introduced to the solid lines, but we found that their contribution is insignificant, at least for the lowest order 'self energy'. Additionally, the resulting equations are more cumbersome than RPA. Note that Dyson equation is usually derived for quantum systems and it is not clear if it holds for our simplified classical perturbation theory.
To make sure that equation (33) is correct, in Hartree-Fock diagrams we replaced single dashed line with double solid line, and summed the resulting infinite series of diagrams. After some combinatorial calculations, we obtained the same expression as in the exponent of equation (33).

Two-dimensional potential
First, we demonstrate our theory with a two-dimensional potential of the form . 34 x y x y 0.05 1.6 5 4.8

2
We deliberately chose specific parameters to make the recrossing effects pronounced. Our test consists of two parts: first we test the equation of motion of a TS trajectory, and then we test the analytic expression for κ. A contour plot of the potential energy surface is shown in figure 3(a). From this plot it is apparent why a large fraction of the trajectories recross the TS. The ridge (black line) is tilted with respect to the TS (blue line). Even after crossing the TS, many trajectories still needs to go uphill to cross the ridge. Thus, there are significant corrections to the x = 0 TS. These corrections, quantified by κ in equation (19), are based upon our analytic determination of which trajectories are reactive. From equation ( We can test this condition by numerically solving the equation of motions. Specifically, we place an atom at the origin with a set of initial velocities. By varying the velocity, we can find the conditions for which the particle oscillates around the TS for a long time. Figure 3(b) compares our numerical solutions with equation (19) for y=0. We see good agreement at the low velocities, which are most important, and only small deviations from the theoretical values at higher velocities. Next we take an ensemble average to compute κ. We sample initial positions at the TS using a rejection sampling method. Once we have random initial positions and velocities, we determine whether trajectories cross or recross using two approaches. The first is the theoretical prediction based on equation (19). This method is partially numerical (sampling initial positions) and partially analytic (using equation (19). The second method is based purely upon simulation: we run trajectories using the Verlet algorithm [24] until the atom moves far away from the TS. If a trajectory is successful, we add its initial velocity to our cumulative flux. At the same time, we compute the TST flux by summing all the positive velocities. The ratio of the true flux and the TST flux is the transmission coefficient, κ. For the calculation of κ we considered several thousands of trajectories.
The results are compared in figure 3(c). We can see that the RPA approximation captures the reduction in κ but there are noticeable differences as compared with the numerical simulation. The non-monotonic shape of the numerically computed k ( ) T function, especially around 400 K, is associated with an abrupt change in the shape of the TS away from the saddle (not visible on the figure). Thus, the harmonic approximation of equation (20) introduces systematic errors. In addition, we also see a discrepancy between the purely numerical simulation and the 'theory + numerics' results, especially at high temperature. This error is caused by limitations of equation (19) which describes the critical conditions for the TS trajectory based upon the interaction Hamiltonian of equation (8)), where we ignored higher order interactions, such as Q Q Q Q i j k 0 . Therefore, we systematically underestimate the recrossing effects. For higher order Hamiltonians, we would also have higher order terms of X in equation (19). When we perform the ensemble average, terms analogous to the ) k T B 2 . A numerical determination of these higher order interaction terms is computationally very expensive. Fortunately, as we will see in the next section, for more realistic (and smooth) systems, errors due to higher order coupling are smaller.

Al surface
The two-dimensional example was considered to illustrate our method for the simplest possible case and to demonstrate the theory with a geometric interpretation. In this second example, which is more realistic, we study different diffusion mechanisms of Al adatom on the Al(100) surface. The interaction between Al atoms is modeled with an embedded atom model [25]. The diffusion mechanism we consider are the following: single atom hop, two-atom, and four-atom exchange. In the two atom exchange, an adatom is pushed into the surface, takes the place of one of the surface atom, and this surface atom pops up on top of the surface in a neighboring site. The four-atom exchange mechanism is similar (figure 4) but with four atoms participating.
The four-atom exchange has a significant dynamical correction factor so we chose to study it in detail (see figure 4) considering different levels of our theory. First, a purely numerical approach was used to sample initial positions and velocities at the TS, as in the previous section. From each initial condition, we ran Newtonian dynamics to determine whether the crossing leads to a reactive trajectory or to a recrossing of the TS. Initial positions were sampled using independent Metropolis-Hastings calculations. As a jumping function, we use (we refer to this method as MHS). The second method is hybrid of our analytical and numerical simulation. We generate initial positions and velocities using the same method, however, we determine the crossing possibilities analytically (we refer to this method as MHT). MHT is a numerical integration of equation (20) without assuming the harmonic approximation. The third method is the same as MHT except that we sample initial positions using the harmonic approximation (here called the HT method). This is equivalent to solving equation (26) numerically. What makes the last method particularly useful is that we do not need to compute energies of the actual system, as it is needed for first two methods. Therefore, the HT approximation to κ can be done much faster than first two methods.
In figure 4(c) we compare these three methods, including those described with Hartree-Fock and RPA. RPA shows remarkable agreement with MHT and HT. There is almost no difference between HT and MHT, which means that there is no need to use Monte-Carlo sampling for the initial positions; the harmonic approximation is sufficient. Most importantly, all three methods agree reasonably well with the direct simulation (MHS).
Results for all three Al diffusion mechanisms at T = 500K are summarized in table 1. RPA works very well for all cases. Note that the Hartree term is zero within numerical error. The reason for this is due to the symmetry of the TS with respect to the reactant and product states. For non-symmetric reactions we may need to introduce higher order diagrams or screening to x x á ñ i j 0 correlation line. HT is also very reliably, with errors between 10 and 20%.

Discussion and Conclusions
In this paper we have introduced a correction to TST based upon stochastic interactions in the vicinity of the TS. For adatom diffusion on an Al(100) surface we showed that TST can overestimate the true rate by as much as 50% making it necessary to compute the dynamical correction factor if an accurate result is required. Most importantly, our approach has the same order of magnitude computational cost as HTST and is almost as accurate as a direct numerical simulation. A summary of our algorithm for computing the transmission coefficient is the following: • Compute the Hessian at the TS.
• Determine the derivative of the Hessian matrix along the negative mode at the TS using, for example, a central difference approximation.
• To compute transmission coefficient, use either the RPA approximation or the harmonic sampling method in combination with equation (19).
The most costly part of these computations are the first two steps, but they can be done for modestly sized systems even when the forces are calculated at the quantum level. The rest of the computation can be done with negligible computational work, in comparison. One of the strengths of our analytic approach is that it provides an extension of TST to semi-classical theory. We have found a mathematical connection between dynamics at the TS and quantum many body theory. If we replace ξ, which are position and momentum variables, with the corresponding quantum operators in equation (26), we obtain a semi-classical theory. It will not be fully quantum theory because the trajectories that we obtained are within the classical limit.
In addition to avoiding the high computation costs of running trajectories to calculate the transmission coefficient, our theory provides new insight into recrossing phenomenon. In equations (30) and (31) both the Hartree and Fock terms are positive and they increase with the number of vibrational modes. However, the derivative of the Hessian is small, and the contribution from most of the vibrational modes is small. For a qualitative estimate we can assume that the motions corresponding to significant displacement contribute most significantly. For the hopping mechanism on the Al surface there is only one atom that is significantly displaced, the other atoms have only small displacements in response to the diffusion mechanism. For the two-atom exchange only two atoms move and so on. This picture is in agreement what we obtain from direct simulation table 1. Therefore, we expect significant recrossing effect when there is a collective motion of many atoms in the reaction mechanism.
Our results are in qualitative agreement with Kramers [40] and Grote-Hynes (GH) [30] theories. According to Kramers and GH theories, the transmission coefficient decreases as the imaginary frequency (w 0 ) decreases and the friction increases. To see that we have obtained a consistent result, note that both a H and a F are proportional to the square of coupling constant Λ. This coupling constant quantifies the interaction strength between the reaction coordinate and the remaining coordinates. In comparison to the GLE, Zwanzig showed that the friction can be expressed using similar coupling constants within a model potential that describes a Brownian particle in a thermal bath [19,32]. For our system too, we can take L 2 as a measure of the friction. Then, according to equation (33), as friction increases, the transmission coefficient decreases, as in Kramers and GH theories. We can also see from equations (15)-(18) that the recrossing effects will be stronger for smaller imaginary frequencies at the TS because there is a higher power of w 0 in the denominators of the A, B, and C matrices. The dependence upon the imaginary frequency at the TS again agrees with Kramers and GH theories. An intuitive explanation why smaller frequencies give rise to a lower transmission coefficient is that for slower motion, there will be more random collisions before crossing the TS and therefore higher chance of recrossing. As we know, the random force and the friction are connected to each other through the fluctuation-dissipation theorem [19]. Whether we choose to describe recrossing phenomena in terms of friction or random collisions, the underlying physics is the same.
This work aims to accomplish the same task as VTST without the need to compute the partition function at points along the reaction coordinate. While this is a significant strength, one shortcoming of our theory is that it cannot be used for reactions without a well-defined TS. In some surface desorption reactions there is no saddle point to define a TS, but rather a monotonous increase towards an energy plateau. Further research is needed to generalize the theory to such cases. When there is a well-defined saddle point, however, the analytic dynamical correction factor presented here is a computationally efficient way of improving the accuracy of HTST. Appendix. Derivation of the equations of motion In this appendix we show how we arrived at equation (12), which can be written as Suppose we know Q(t). We can treat Q i (t) terms as an external force We can write the equation of motion in matrix form The general solution of equation (A.4) is We use the following property of matrix A to compute the exponents w = ( ) A I, A . 9 From the initial condition