Reaction rate theory has been at the center of physical chemistry for well over one hundred years. The evolution of the theory is not only of historical interest. Reliable and accurate computation of reaction rates remains a challenge to this very day, especially in view of the development of quantum chemistry methods, which predict the relevant force fields. It is still not possible to compute the numerically exact rate on the fly when the system has more than at most a few dozen anharmonic degrees of freedom, so one must consider various approximate methods, not only from the practical point of view of constructing numerical algorithms but also on conceptual and formal levels. In this Perspective, I present some of the recent analytical results concerning leading order terms in an ℏ2m series expansion of the exact rate and their implications on various approximate theories. A second aspect has to do with the crossover temperature between tunneling and thermal activation. Using a uniform semiclassical transmission probability rather than the “primitive” semiclassical theory leads to the conclusion that there is no divergence problem associated with a “crossover temperature.” If one defines a semiclassical crossover temperature as the point at which the tunneling energy of the instanton equals the barrier height, then it is a factor of two higher than its previous estimate based on the “primitive” semiclassical approximation. In the low temperature tunneling regime, the uniform semiclassical theory as well as the “primitive” semiclassical theory were based on the classical Euclidean action of a periodic orbit on the inverted potential. The uniform semiclassical theory wrongly predicts that the “half-point,” which is the energy at which the transmission probability equals 1/2, for any barrier potential, is always the barrier energy. We describe here how augmenting the Euclidean action with constant terms of order ℏ2 can significantly improve the accuracy of the semiclassical theory and correct this deficiency. This also leads to a deep connection with and improvement of vibrational perturbation theory. The uniform semiclassical theory also enables an extension of the quantum version of Kramers’ turnover theory to temperatures below the “crossover temperature.” The implications of these recent advances on various approximate methods used to date are discussed at length, leading to the conclusion that reaction rate theory will continue to challenge us both on conceptual and practical levels for years to come.

In 1889, Arrhenius published his seminal paper on chemical reaction rates entitled “On the reaction rate during the inversion of cane sugar by acids.”1 As detailed in a review article by Laidler,2 Arrhenius was not the first to write down the rate expression in terms of an exponential inverse relation on the temperature T,
rateAexpEkBT.
(1.1)
The innovation in his paper was the empirical understanding that the energy E in the exponent signifies the energy of an intermediate state, which must be crossed before reaction occurs. This was a seminal step in the research for a microscopic first-principles understanding of the exponential form and the prefactor A.

Tolman, in the closing remarks of his book published in 1927,3 writes as follows: “The problem of reaction velocities is probably nearer to the heart of most chemists than anything else in their whole range of activity. Rates of reaction are the factors that determine yields, and costs, and possibilities, and their theory must eventually succumb to scientific treatment…. Whatever may be, the way is long and the traveler must not tarry.” The history of the development of rate theory in the 20th century may be found in many textbooks and review articles, a definitive one is that of Hänggi, Talkner, and Borkovec.4 As described there, one of the landmarks in the development of rate theory was Farkas’ 1927 paper5 in which he put forth that the basis for a theoretical calculation of rates is the assessment of the ratio of flux of particles going from reactants to products to the density of reactants. This laid the ground work for transition state theory (TST) as derived by Wigner6–8 and postulated by Eyring9,10 and enabled the numerical computation of reaction rates. The relationship between the theories developed separately by Eyring and Wigner is described in Ref. 11.

Less than a decade later, Kramers laid the groundwork for a stochastic theory of rates, in which he considered the interaction of a particle with a medium and the resulting effect of the medium on the rate of reaction.12 Kramers considered two central effects. One is rate determining when the interaction of the system with the surrounding medium is weak, known as the energy diffusion limited regime. The bottleneck is providing sufficient energy for the system to overcome the barrier. The second, known as the spatial diffusion limited regime, is rate determining when the interaction becomes strong and prevents the system from crossing smoothly through the barrier.

One of the challenging questions that plagued rate theory for many years was the accurate interpretation of the exponent appearing in the Arrhenius expression (see especially the discussion on p. 25 of Ref. 13). It was only in 1978 that Chandler established the fact that when considering thermal rates, it is a free energy, not the bare potential energy of the lowest energy barrier separating reactants and products.14 

During the first half of the 20th century, it also became clear that quantum effects cannot be ignored. Roughly, they affect the rate in three ways. One is due to the quantization of the energy levels of the reacting system. Eyring dealt with this difficulty by recasting the rate expression in terms of partition functions so that one could employ quantum partition functions. A second effect has to do with tunneling, a mechanism which remains intriguing to this very day, as may be seen from some of the books published on the topic.15–22 A third effect is that of quantum resonances;23,24 however, typically when considering thermal rates, these are important only when the temperature is very low.25 

How then does one compute the thermal rate? The flux over population expression of Farkas is readily estimated using classical mechanics. One chooses an initial phase space point from the thermal distribution of reactants and propagates a trajectory using classical mechanics, and if it ends up in the product region, it contributes positively to the total reactive flux.26 The classical algorithm was formalized almost simultaneously by Miller27 and McLafferty and Pechukas28 into a quantum expression whereby the flux was considered as the quantum mechanical trace operation of a product of three operators—the thermal operator expβĤ, with β = 1/(kBT) and Ĥ being the Hamiltonian of the system; the quantum flux operator; and a projection operator, which takes into account only those states that are reactive, that is, they lead from reactants to products. Subsequently, Miller, Schwartz, and Tromp, in a now classical paper,29 reformulated the expression by considering the rate as a time integral of the thermal flux autocorrelation function or equivalently the infinite time limit of a flux-side correlation function. This sets the groundwork for the numerically exact computation of thermal reaction rates in complex systems.30 However, the practical challenge is not trivial. To obtain the rate, it is necessary to propagate to infinite time, and doing this quantum mechanically in large systems is even today not feasible.

This is perhaps one of the reasons that Miller developed a semiclassical rate theory, which could, in principle, account for both quantization and tunneling albeit approximately.31,32 A central advantage of a semiclassical theory is that it utilizes classical trajectories, which with the modern day quantum chemistry technology may be computed on the fly. Trouble is that semiclassics has phases so that converging semiclassical rate theory using Monte Carlo methods is also much too expensive.33 All this leads to the happy situation that even today the computation of chemical reaction rates, even when considering motion on only one Born–Oppenheimer potential energy surface, is still an open problem. A central part of this Perspective is to try and understand some of what has been done and consider some aspects of present day theory, which can be improved upon.

One should distinguish between methods that lead to computational algorithms and analytic results. One very important model, which was solved analytically by Bell,34 is that of the thermal transmission through a parabolic barrier. Wigner, in his 1932 paper on tunneling rates,7 derived only the leading order correction term of order ℏ2 to the classical result; however, he already then recognized that even in one dimension, a parabolic barrier is not the whole story. By expanding the potential to fourth order about the barrier top, he derived an analytic expression, correct up to order ℏ2 for the thermal rate of transmission through a potential that is symmetric about the barrier top. It is only very recently that this expression has been generalized to include asymmetric barrier crossings35 and provide the next term of the order of ℏ4.36 Why is this at all of interest? Let us remember that there are many useful approximate ways to compute the thermal rate, but which is the best? Analytic results help in distinguishing between different approximations. We shall also show below how the analytic results may be used to improve upon existing semiclassical approximations.

Related to this is the vibrational perturbation theory (VPT) of rates as developed by Miller et al. in 199032 and further expanded by Hernandez and Miller in 1993.37 Vibrational perturbation theory was developed in the middle of the previous century for spectroscopy.38,39 Using second order quantum perturbation theory, explicit expressions were derived for the quantized energy levels of bound anharmonic systems. Miller et al.32 extended these results to reactions where at the lowest energy saddle point the reaction coordinate is characterized by an imaginary frequency. In this theory, one expands the potential around the barrier top, including up to fourth order spatial derivatives in all degrees of freedom. Then, by assuming a quadratic energy action relation, one obtains an explicit expression for the tunneling action as a function of the energy. This is then used within the context of a multidimensional semiclassical approximation for the tunneling probability.31 Barker and co-workers40,41 further developed this approach, including in it the anharmonic correction for the ground state energy, to derive what is known as VPT2 theory. Since it needs as input only the derivatives of the potential at the barrier, it is more cost-effective to implement than running trajectories on the fly.42,43 However, even in one dimension, the theory fails for asymmetric potentials.44 Can one improve upon VPT2 theory without paying a heavy computational price?45 

Many approximate theories have been developed in parallel to the development of numerically exact quantum algorithms.46,47 The central demands from a practical algorithm are that it would be amenable to on-the-fly computations, that it would account in a reliable manner for the effect of quantization of molecular levels, and that it includes quantum phenomena, such as tunneling. In this context, one should mention the classical Wigner dynamics method,48–51 initial value semiclassical methods,52–57 Centroid Molecular Dynamics (CMD),58,59 Ring Polymer Molecular Dynamics (RPMD),60–62 semiclassical instanton theories,31,63–71 coupled coherent states methods,72 and quantum instanton approaches.73–76 As this is a perspective article, we will not review here in detail each of these methods, but will mention them wherever relevant, attempting to stress their advantages and shortcomings, always with the goal of pointing out future ways of improvement.

The formalism considered thus far is, in principle, applicable to any system. However, reactions in condensed phases present a special challenge. Even within the realm of classical mechanics, it is not yet possible to propagate on-the-fly the many thousands of degrees of freedom involved in a reaction occurring in a condensed phase.77,78 Even if it were possible, the results of such ab initio molecular dynamics computations would be rather opaque. A simplifying but realistic model is very helpful for making sense of measured phenomena and guiding our understanding of the mechanisms. Arguably, the best model for such purposes goes under the terminology of “dissipative systems” in which the “system degrees of freedom” are coupled bilinearly to a bath of harmonic oscillators.79 In classical mechanics, this is equivalent to a generalized Langevin set of equations, in which the system degrees of freedom are affected by the frictional and random forces induced by a bath.80 

The dynamics of dissipative systems has a long history, initiated by the seminal paper of Kramers. For many years, Kramers’ theory and transition state theory were considered to be separate entities. It was only in the 1980s that it was established that in the spatial diffusion limited regime, Kramers’ rate expression and its generalization to include memory friction effects81,82 is identical to transition state theory83 since the Langevin equation underlying Kramers’ theory may be derived from a Hamiltonian. This also led to the observation that Wolynes’ quantum analytic rate expression for a dissipative parabolic barrier84 may be derived from a quantum version of transition state theory.85 At the same time, one of the challenges presented by Kramers, known as the Kramers turnover problem, was to derive an analytic expression for the rate valid for all friction values. This challenge was answered by Mel’nikov and Meshkov86 and formalized into what is known today as Pollak, Grabert, Hänggi (PGH) theory.87 Yet, there remained many open questions, foremost among them was the formulation of a quantum version of PGH theory,88,89 which was not limited to parabolic barrier based quantum dynamics.

The turnover theory is not just a purist question of finding an analytic solution to a challenging problem. Kramers’ turnover has been observed experimentally using levitated nanoparticles.90 There have been reports of its observation for reactions in liquids.91–95 It is especially applicable to surface diffusion, where even hopping distributions have been measured96–98 and interpreted theoretically using the turnover formalism.99 

This Perspective is organized as follows. To set the stage, we briefly review in Sec. II the correlation function formalism29 for computing reaction rates and the various approximate methodologies developed in recent years for implementing the formalism to real systems. Section III presents recent analytic expansions of the transmission coefficients in terms of ℏ2n35,36 and especially their implications for the vibrational perturbation theory approach to computing reaction rates.100 This also serves as an introductory to Sec. IV where tunneling theories are considered, paying special attention to instantons and their role in thermal and microcanonical contexts. Recent developments in Kramers’ turnover theory101–103 are considered in Sec. V. We end with Sec. VI discussing the open problems, which we foresee for the coming years.

The fundamental expression underlying the computation of thermal rates is the flux to population ratio. Without loss of generality, we assume N + 1 degrees of freedom, q denoting a reaction coordinate and x denoting all the others. The potential energy surface is characterized by a saddle point located at q,x̱. The population of reactants is assumed to be thermal and given by the partition function of reactants,
Za=TrexpβĤθqq,
(2.1)
where θx is the unit step function, β = 1/(kBT) is the inverse temperature, and Ĥ is the Hamiltonian operator. The flux operator measures the number of particles per unit time crossing a point along the reaction coordinate. For simplicity, this point is taken to be the barrier location q so that the operator is
F̂=12Mp̂qδqq̂+δqq̂p̂q,
(2.2)
where “hats” denote operators, M is the mass associated with the “reaction coordinate,” and p̂q denotes the momentum operator conjugate to q̂. A central part of the formulation of Miller et al.29 was to realize that the projection operator that chooses only the reactive contribution to the flux has the form
P̂=limtexpiĤtθq̂qexpiĤt.
(2.3)
It is the time evolved spatial step function that guarantees that only contributions that at long time come from evolution from the reactant to the product side contribute to the rate. One may equally well replace, and this was the original formulation,28,29 the spatial step function with a time evolved momentum step function, but the spatial form is the more computationally convenient one, even though it leads to a vanishing flux in the zero time limit,104 while the momentum version does not.
The rate Γ is given by the ratio of the reactive flux to the population of reactants, and so it takes the form
ZaΓ=TrexpβĤF̂P̂=TrexpβĤ2F̂expβĤ2P̂,
(2.4)
where the second equality results from the fact that the projection operator commutes with the Hamiltonian. One may use this to rewrite the expression in a Kubo form,105–108 which has some practical advantages. The analogous classical expression for the rate is
Za,clΓcl=limtdpqdq dp̱xdx̱2πN+1expβHclpqMδqqθqt,
(2.5)
where Hcl denotes the classical Hamiltonian and qt denotes the time evolved value of the reaction coordinate of a trajectory initiated at the phase space point pq,q;p̱x,x̱. Za,cl is the classical analog of the reactant partition function. A central quantity that will be of interest is the quantum transmission coefficient, defined as
κβ=ZaΓZa,clΓcl,
(2.6)
expressing the ratio of the quantum reactive flux to the classical. Both quantum and classical reactive fluxes are independent of the location of the point at which they are measured;26 this is an expression of the conservation of the number of particles.
The flux side expression may be recast formally using coherent states109 into an illuminating form. In one dimension, a coherent state is just a Gaussian,
x|gp,q=γπ1/4expγ2qx2+ipqx.
(2.7)
Since the reactive flux is independent of the precise location of the dividing surface, one may use a “smeared flux operator.”28 With a Gaussian smearing function, one can then show that (written here in a one-dimensional form)
ZaΓ=0dxdpdq2πpMδqq×xexpβ2+itĤgp,q2.
(2.8)
This expression has a few advantages, as the averaging (apart from the velocity factor) is over a positive definite quantity, and it has a classical flavor to it. It would be natural to employ the coupled coherent states propagation method of Shalashilin and Child72 for the complex time evolution. It also expresses the fact that, in principle, only one time propagation is needed instead of both forward and backward time propagations as implied by Eq. (2.4). This reduces the numerically exact effort and can be used to simplify semiclassical approximations. Its full use remains to be explored.

The central objective of transition state theory (TST) is to provide a thermodynamical approximation to the reactive flux, which does not demand any real time propagation. In classical mechanics, this is straightforward, and one introduces the concept of a dividing surface, defined so that all trajectories going from reactants to products must cross it,26 and approximates the projection operator to be θpq, that is, one chooses only those trajectories moving toward products as they cross the dividing surface. It is well known that this bounds the reactive flux from above, leading to variational transition state theory, whereby one searches for the surface that divides between reactants and products, which gives the lowest unidirectional flux.110,111 Significant progress has been made in creating algorithms for computation of optimal dividing surfaces.112–115 However, Wigner already noted that a quantum analog to classical transition state theory is a difficult proposition, as it demands the simultaneous knowledge of the location of the particle and the sign of its momentum, and this is disallowed by the uncertainty principle.8 This does not mean that one cannot formulate an approximation that is thermodynamical in nature, which will give reasonable results. However, any attempt at bounding the reactive flux from above has failed in the sense that rigorous upper bounds are not very close to the exact result even for a parabolic barrier potential.28,116,117 For example, the coherent state rate expression of Eq. (2.8) will give a strict upper bound to the reactive flux by limiting the momentum integration to 0,; however, the resulting upper bound, in one dimension, is identical to the upper bound of McLafferty and Pechukas,28 and as already mentioned, it is not very accurate.

One way of implementing a quantum thermodynamic strategy is by first recasting the reactive flux expression as a phase space trace operation in Wigner phase space.118,119 The Wigner representation of an operator Ô in phase space is
Opq,q;p̱x,x̱=dξqdy̱expipqξq+ip̱xy̱×qξq2,x̱y̱2Ôq+ξq2,x̱+y̱2.
(2.9)
One readily finds that the phase space representation of the flux operator is the same as the classical. One then uses the exact quantum representation of the thermal and flux operators but replaces the quantum projection operator with its transition state theory form θpq. As discussed in Sec. III B, this involves two approximations. One is ignoring the fact that the trace of three operators is not simply the phase space trace of the product of the three Wigner representations of the operators. The second is the TST replacement of the exact projection operator.
These approximations notwithstanding this approach lies at the heart of the semiclassical approximation to the reactive flux. Miller has shown31 that this leads (written here for the sake of brevity in one dimension) to Kemble’s expression120 for the energy dependent transmission probability,
TuscE=11+expSE,
(2.10)
where (in one dimension)
SE=2qEq+Edq2MVqE
(2.11)
is the Euclidean action on the upside down potential energy surface with q±E denoting the two turning points defined by Vq=E. The subscript usc in which the letter u has been added to the transmission probability in Eq. (2.10) is to remind us that this expression is the uniform generalization of the primitive semiclassical approximation,
TscE=expSE.
(2.12)
It is interesting to note that although Miller’s derivation at the outset was based on a quantum transition state theory approximation of choosing only positive momenta, the microcanonical result in one dimension [Eq. (2.10)] was originally derived by Kemble120 and formalized by Fröman and Fröman,121 without resort to any transition state theory assumption. It was derived semiclassically by imposing suitable uniform continuity conditions at the classical turning points.

Miller’s derivation was not limited to one dimension. It related tunneling even in many-dimensional systems to a classical object—a periodic orbit moving on the inverted potential. The thermal transmission coefficient is obtained by averaging over the energy dependent transmission coefficient with the canonical weight expβE. Miller used the semiclassical transmission probability TscE [Eq. (2.12)] in conjunction with a steepest descent estimate of the integral over the energy, to derive a general semiclassical theory for thermal rates, centered about the “instanton” periodic orbit moving on the inverted potential. A few years later, Callan and Coleman122 rederived Miller’s thermal expression, by semiclassical estimation of the imaginary part of the free energy, a method that, in principle, is exact20,64,79,123 at sufficiently low temperature (contrary to a dissenting claim124). Only recently, Miller’s derivation of the multidimensional energy dependent transmission coefficient has been formalized on the basis of the microcanonical flux autocorrelation function.125 

A different class of approximations is based on treating the thermodynamic part numerically exactly and only replacing the long time quantum limit of the projection operator with a classical long time limit. This may be implemented in different ways. Cao and Voth suggested replacing the quantum density operator with a centroid density and then propagating classical trajectories on the resulting centroid density effective potential.58,59,126,127 This approximation is known as centroid molecular dynamics (CMD). A different suggestion was to use the exact phase space representation of the symmetrized thermal flux operator but treat the time evolution of the projection operator classically; this is known as mixed quantum classical theory,49,108,128–131 evolving from an earlier suggestion made by Heller.48 It also goes under the nomenclature of linearized semiclassical theory.52 It is exact for short times.132 A third suggestion was to replace the Hamiltonian appearing in the rate expression by its ring polymer form and then treating the flux side correlation function as a classical object, albeit in the enlarged ring polymer phase space. This approximation is known as ring polymer molecular dynamics (RPMD).60–62,133,134 These approximate methods have been reviewed elsewhere, so here, we do not go into their details. Each of these has advantages and disadvantages, some of which have become obvious only lately, and they will be considered in Secs. III, IV, and VI.

There are only a few one-dimensional models for which one may derive the exact expression for the transmission coefficient. Arguably, the most influential one is the parabolic barrier potential,
Vq=VMω2q22.
(3.1)
In one dimension, the exact expression for the one-dimensional thermal transmission coefficient is
κβ=βexpβV0dEexpβETE.
(3.2)
For the parabolic barrier, the exact energy dependent transmission coefficient is well known,120 
TpbE=1+exp2πωVE1.
(3.3)
The energy interval for the parabolic barrier potential is ,, while in reality, there is always a lower bound to the energy; without loss of generality, it is taken to be E = 0 for reactants.
Inserting the parabolic barrier transmission probability into the expression for the thermal rate and changing variables from E to x=βEV give
κpbββVdxexpx1+exp2πxu1,
(3.4)
where we used the notation
u=βω.
(3.5)
In many cases, the reduced barrier height βV ≫ 1, so one may replace the lower limit −βV with −∞ to obtain the celebrated parabolic barrier exact expression for the thermal transmission factor,34 
κpbβ=u2sinu2.
(3.6)
This result diverges at the so-called (inverse) crossover temperature,
βc=2πω.
(3.7)

The divergence is obviated by noting that if u > 2π, the exponent exp2πu1x will grow without bounds as x → −∞. Thus, for “high” temperatures defined by u < 2π, the parabolic barrier approximation is useful, and in this temperature region, energies above the barrier contribute mostly to the rate. At lower temperature, the parabolic barrier approximation fails, and one must consider anharmonic potentials for which the rate does not diverge, and in this regime, the main contribution to the thermal rate comes from energies that are below the barrier, thus the nomenclature “crossover temperature;” it signifies the crossover from a tunneling mechanism to a classical-like above barrier transmission mechanism.

Practically, this limitation is quite meaningful in the context of chemical reactions. A barrier frequency of n · 100 cm−1 with n = 1, 2, 3, … implies a corresponding “crossover” temperature of Tnn·24 K, implying that at room temperature many reactions, especially those that are characterized by heavy atom transfer, will not exhibit deep tunneling effects. It is then not surprising that in the “heavy atom community,” an accepted way of accounting for quantum tunneling effects is to use the parabolic barrier transmission coefficient [Eq. (3.6)] or its high temperature limit of 1 + u2/24.135,136

Analytic expressions for the exact quantum energy dependent transmission coefficient and the instanton action in one dimension exist for the Eckart,137 inverse Morse,138 and square barrier potentials.139 However, for all these, the thermal rate cannot be obtained analytically exactly; at best, one may derive an analytic steepest descent estimate based on the known energy dependent result.

Even derivation of analytic expressions in the high temperature limit is challenging. Wigner in his 1932 paper7 showed that the leading order correction term to the transmission coefficient is not merely the parabolic barrier result of 1 + u2/24, but there is an additional term, which for a symmetric potential is proportional to the fourth derivative of the potential (V4) at the barrier. The exact ℏ2 correction term valid also for asymmetric potentials
κβ=1+u2241+14V4βV22V323βV23
(3.8)
has been derived only recently35 (here, Vj denotes the j-th derivative of the potential at the barrier top). In most chemically relevant cases that we are aware of, the fourth derivative is positive, but for example, in the inverse Morse potential, it is negative, implying that at high temperature, the quantum correction can be dominated by above barrier quantum reflection, leading to a thermal transmission coefficient, which is less than unity. Approximate theories, such as CMD, RPMD, and mixed quantum classical Wigner dynamics, do not account for this. Recent (controversial140–142) attempts at formalizing RPMD as “the” quantum transition state theory143–145 also seem to be insufficient, as they do not give the correct leading order correction term. At best, they give the correct parabolic barrier leading order correction. Although heavy atom transfer reactions would seem to be classical in nature, justifying the use of these three approximations, it turns out that in reality they may be insufficient. For example, for the Eckart barrier, the parabolic barrier approximation underestimates the exact rate.35 The only redeeming aspect is that the error involved in the approximation will typically not be more than 50%, and this is within the typical accuracy of rate constant measurements.

The derivation of the high temperature limit in one dimension turns out to be straightforward though quite involved. As already noted, the quantum reactive flux expression is a trace of a product of three operators. In one dimension, their phase space representation in the classical limit is well known. The thermal density is expβHcl, the flux is pqMδqq, and the projection operator involves choosing only those values of the momentum, which lead to energies that are greater or equal to the barrier height. For reactions that proceed from left to right, the sign of the momentum is positive so that the classical projection operator at any point in phase space in one dimension is49  θp2MVVq, where the positive (negative) sign is to be taken for qq (qq).

To obtain the leading order quantum contribution, one uses the leading order quantum correction term for the density, which already appears in Wigner’s 1932 paper.118 Then, one must obtain the leading order correction term to the projection operator, and this may be found in Ref. 49. Since the trace involves three operators, one must obtain the leading order correction to the phase space distribution of a product of them. In general, when considering operators  and B̂ and their exact Wigner representations Ap,q,Bp,q as defined in Eq. (2.9), the phase space representation of their product Ĉ=ÂB̂ is119 
Cp,q=Ap,qexpi2Λ̂Bp,q,
(3.9)
where the Janus operator Λ̂ is defined as
Λ̂=qppq.
(3.10)
The arrows define the operator on which to apply the corresponding derivative. To obtain the second order term in ℏ, one expands expi2Λ̂1+i2Λ̂28Λ̂2+. For our purpose, the linear term in ℏ will vanish since the reactive flux is real. Putting all this together, with some algebra, gives the leading order correction [Eq. (3.8)].

One need not stop here; the explicit expression for the next term that goes as ℏ4 was derived in Ref. 36. The analytic result is not trivial and includes up to the eighth derivative of the potential at the barrier. At present, this makes it prohibitive to compute on the fly unless one has an analytic expression for the potential energy surface.

Is all this “important” or is it just an algebraic exercise? As already mentioned, it offers a test of approximations, but as we shall see below, and perhaps this is not less interesting, it also serves as a guide for improving the semiclassical approximation for the energy dependent transmission coefficient. A missing ingredient in all of this is the generalization to many dimensional systems. Apart from the complexity of the resulting expressions, there is a much more fundamental difficulty. Unless the potential is separable, the classical limit of the projection operator is not known analytically. Due to the interactions between modes, trajectories may recross the surface at which one measures the reactive flux. Only by propagating forward and backward in time, can one know whether they contribute to the reactive flux. Does this imply that it is impossible to obtain even the leading order term in ℏ2 in the multidimensional case? Not necessarily, for example, one may use a perturbation theory for the classical trajectories in which the anharmonic part of the potential serves as the “small parameter.”49 However, this remains to be worked out; if you like analytic manipulations, go ahead and do it. It is “important.”

To present the central ideas underlying the theory, it is useful to consider it in one dimension. For the sake of completeness, the multidimensional result will be presented without going into a detailed derivation. The “standard” semiclassical transmission coefficient has been given in Eq. (2.12). Especially in the deep tunneling region, the action SE is large compared to ℏ, but its magnitude reduces as the energy increases so that the transmission coefficient increases with the increasing energy. On the other hand, especially at low temperature, the exponential term expβE decreases exponentially with the increasing energy. The product of the two terms is then rather narrow in an energy bell shaped function, making it an ideal candidate for estimation via steepest descents. Since the action is on the inverted potential, the period of the instanton is given by
τE=dSEdE.
(4.1)
The stationary point of the exponent βE+SE/ is thus given by
β=τEβ.
(4.2)
This identifies the thermal instanton as the periodic orbit on the inverted potential with period ℏβ and energy Eβ. Expanding the exponent to second order in the energy and performing the Gaussian integral by letting the lower limit of the energy go to −∞, one readily finds the result
κscββexpβV2πS2EβexpβEβSEβ,
(4.3)
where S2Eβ is the second derivative of the instanton action with respect to the energy, at the steepest descent energy Eβ. It is not necessary to extend the integration to −∞; keeping the finite lower limit just adds an error function term to the estimate.64 

As the temperature increases, the period of the instanton decreases and the energy Eβ increases. At some temperature, the period ℏβ will equal the period of the instanton of the parabolic barrier, which is 2πω. This will occur when EβV. The period of a harmonic oscillator is energy independent so that the second energy derivative of the action S2Eβ, which is the first energy derivative of the period τ(Eβ), will vanish and the transmission coefficient will diverge. As seen, this occurs when ℏβω = 2π, and we again meet the “crossover temperature,” which now gets a new meaning. Below the crossover temperature, there exists an instanton orbit and transmission occurs by tunneling, guided by the instanton. Above the crossover temperature, transmission occurs mainly at above barrier energies. In reality, even at the barrier energy, the second derivative does not necessarily vanish,146,147 but it will be small, leading to an estimate, which is much too large.126 Attempts have been made in the past to overcome this issue,64,65,67,126 but in the following, we will demonstrate that when using the uniform semiclassical theory, the problem naturally disappears.71 

The multidimensional version of the steepest descent estimate of the rate for a system with N + 1 degrees of freedom is
ZaΓ12πβπj=1N12sinhβωjEβ2πβ22S2Eβ×expβEβSEβ,
(4.4)
where the frequencies ωjEβ,j=1,,N, are the stability frequencies of the instanton orbit at the steepest descent energy Eβ and Za is the partition function of reactants. The steepest descent equation for the energy Eβ has the same form as in the one-dimensional case.

The instanton method for estimating the thermal transmission coefficient has many advantages. Foremost, there is no averaging involved. All that one needs is to find the instanton orbit with the least action whose period is ℏβ. Methods have been developed for this purpose65,70,148 so that this is arguably the cost-effective way of obtaining especially deep tunneling rates in multidimensional systems on the fly.68,149,150 The result also takes into consideration the quantization of all stable modes, through the stability frequencies of the orbit, albeit on a harmonic level. The limitation to harmonic quantization is not too drastic, as tunneling is usually found at low temperatures. We note though that an advantage of methods such as mixed quantum classical Wigner dynamics, CMD, and RPMD is that the quantization is not limited to harmonic approximations.

How accurate is the old version of instanton theory? It has been tested extensively in low-dimensional systems. Its accuracy is typically not worse than 50%, and usually much better, especially in the deep tunneling region. When coming close to the crossover temperature, it fails. Care must be taken when the temperature is too high; however, this is the parameter region where the other methods become rather accurate, and so one may use them instead.

One of the central challenges of all approximate methods is to improve them systematically. As of the writing of this Perspective, a practical way of implementing such a strategy for the CMD and RPMD methods does not exist. It can be implemented, in principle, for the classical Wigner approximation49 but, in practice, is much too expensive. The same goes for semiclassical initial value representation methods, for which one may show that there exists a time dependent perturbation series, which converges to the exact answer;151,152 however, going beyond the leading order term is a heroic effort,33,153 even in small systems. Are there ways of improving the instanton method systematically? One may further expand the action about the steepest descent energy to more than second order and obtain corrections.154 Obtaining numerical derivatives of the action that are higher than the second order is though a demanding task.154 Similarly, an assumption in the derivation of the expression is that the stability frequencies of the instanton do not change rapidly with energy. This assumption is, in principle, not necessary as one may follow the energy dependence of the stability frequencies. This would also give an expansion about the instanton energy and could be included in the steepest descent estimate. The obvious difficulty is that it requires a much larger numerical effort, especially when considering on-the-fly computations in large systems. As far as this author is aware, such an expansion has not been implemented even for relatively small systems where an analytic potential energy surface would be available.

There are other challenges. Especially at low energies, quantum resonances can affect the rate.155 The semiclassical approximation underlying the instanton method only includes the tunneling action; phases that interfere and lead to resonances are not included in it. Another low temperature problem has to do with quantum threshold reflection.156 At threshold, unless some unusual resonance conditions hold,157 the quantum transmission coefficient vanishes; this is known as quantum threshold reflection. The semiclassical transmission probability given in Eq. (2.10) will be exponentially small but finite, implying that the instanton theory cannot account for quantum threshold reflection. This is not a difficulty when considering the decay of a meta-stable state, but it does arise in bimolecular reactions. Perhaps though, the most challenging aspect of the instanton theory is its divergence at the crossover temperature. In Subsection IV B, we will describe a recent advance,71 which answers this challenge.

As in Subsection IV A, the description here will be mainly for one-dimensional systems;71 the multidimensional result158 will be presented without derivation. Many of the defects of the “old” instanton theory are removed with a simple change. One may formally rewrite the thermal transmission probability using Kemble’s expression120 for the energy dependent transmission probability [Eq. (2.10)] rather than the primitive semiclassical form of Eq. (2.12). The uniform semiclassical thermal transmission coefficient is then given by
κuscβ=βexpβV0dEexpΦE
(4.5)
with
ΦE=βE+ln1+expSE.
(4.6)
The energy integration in Eq. (4.5) may be estimated using steepest descents.
This is not a minor change. Foremost, the steepest descent energy is now determined somewhat differently. Setting the energy derivative ΦE=0 leads to the steepest descent equation,
β=11+expSEβτEβ.
(4.7)
At low temperature, due to the fact that the instanton action S(Eβ) becomes large relative to ℏ, this condition reduces to the same as in the “old” theory [Eq. (4.2)]. The central difference may be seen when the energy of the instanton equals the barrier energy. At this energy, the instanton reduces to a point (the barrier top) and its action vanishes; the period of the instanton will be the parabolic barrier period 2π/ω. However, due to the denominator, which is now 2, this will occur when ℏβω = π rather than 2π. If one considers the crossover temperature as that temperature at which the instanton energy equals the barrier energy, then here we find that this crossover temperature is a factor of 2 higher than in the old theory.
Not less important is that we will not find any divergence. For the steepest descent estimate, one needs the second derivative of the exponent ΦE at the instanton energy Eβ. This second derivative now has two contributions,
Φ2Eβ=d2ΦEdE2|E=Eβ
(4.8)
=11+expSEβd2SEdE2|E=Eβ+β2expSEβ|E=Eβ.
(4.9)
The first term on the right-hand side depends on the second derivative of the action, and it may vanish at the crossover temperature. However, the second term depends on the period of the action squared, and this does not vanish when Eβ = V.
The steepest descent estimate of the transmission coefficient is
κuscββexpβV2πΦ2EβexpΦEβ.
(4.10)
It is continuous as a function of β and does not suffer from any divergence. It is then of interest to study this expression for the parabolic barrier potential for which the action is
SE=2πωVE.
(4.11)
The solution for the steepest descent energy is
Eβ=V+ω2πln2πu1
(4.12)
with u = ℏβω. When the temperature is sufficiently high, u → 0 and the energy Eβ goes to ∞ as −ln u. At the crossover temperature u = π, the energy Eβ = V, and as u approaches 2π from below, the energy Eβ goes to −∞, and this is what one should expect.
Consider then the steepest descent estimate for the transmission coefficient. The second derivative of Φ(Eβ) for the parabolic barrier is
Φ2Eβ=β22πu1
(4.13)
and
κusc,pbβ2π1u2π12u2π2πuu2π+12.
(4.14)
One notices that in the high temperature limit, the steepest descent estimate is such that
limu0κusc,pbβu0,
(4.15)
and this is clearly incorrect; in the high temperature limit, the transmission coefficient should go to unity.
What went wrong? In the steepest descent approximation, the integrand is a Gaussian in the energy. However, in the high energy limit, we know that the energy dependent transmission coefficient goes to unity so that the integrand should decay exponentially as expβE, and this is also true if the action is different from the parabolic barrier action. It is straightforward to remedy this failure. One should force the Gaussian whose exponent is ΦEβ+Φ2EβEEβ22 to turn into an exponential when the energy is sufficiently high. To do this, one looks for an energy E* such that above this energy, the exponent becomes linear and the general form of the exponent is
ΦE=ΦEβ+Φ2EβEEβ22×θE*E+βE+UβθEE*.
(4.16)
For this to make sense, we want the exponent and its first derivative to be continuous at E = E*, and these two conditions will determine the values of E* and Uβ. It is straightforward to solve the two resulting equations, and one finds
E*=Eβ+βΦ2Eβ
(4.17)
and
Uβ=β2Φ2Eβ1βln1βτEβ.
(4.18)
One notices that the transition to an exponential dependence on the energy occurs at energies that are larger than the steepest descent energy, and this is as it should be. Using Eqs. (4.16)(4.18), the integration over the energy is straightforward and one finds the central result
κuscexpβVEβ1+expSEβπβ22Φ2Eβ×1+erfβ22Φ2Eβ+expβ22Φ2Eβ.
(4.19)

As before, it is of interest to study this result for the parabolic barrier. Now, in the limit that u → 0, the transmission coefficient goes to unity. In Fig. 1, we plot the ratio of this steepest descent estimate to the exact parabolic barrier result [Eq. (3.6)] in the region of 0 ≤ u < 2π. One notes that except when coming close to the “old” crossover temperature (u = 2π), the uniform instanton estimate for the parabolic barrier is within 20% of the exact result. Without the uniform semiclassical theory, such a comparison would not have been possible, indicating the power of the uniform semiclassical instanton theory.

FIG. 1.

The ratio of the steepest descent estimate to the exact result for the thermal transmission coefficient of a parabolic barrier is plotted as a function of the reduced inverse temperature u = ℏβω. The solid line at P = 0.8 indicates the range of temperatures at which the steepest descent estimate is “reasonable.”

FIG. 1.

The ratio of the steepest descent estimate to the exact result for the thermal transmission coefficient of a parabolic barrier is plotted as a function of the reduced inverse temperature u = ℏβω. The solid line at P = 0.8 indicates the range of temperatures at which the steepest descent estimate is “reasonable.”

Close modal
The generalization to the multidimensional case is straightforward,158 and one finds
ZaΓ12πβπj=1N12sinhβωjεβ2×expΦεβ+βVπβ22Φ2εβ×1+erfβ22Φ2εβ+expβ22Φ2εβ
(4.20)
with
ε=Ej=1Nnj+12ωjE
(4.21)
and
Φε=βε+ln1+exp1Sε.
(4.22)
As before, ωjE are the stability frequencies of the multi-dimensional instanton whose energy is E.
The thermal instanton expression, in principle, will not be any better than the numerically integrated energy dependent transmission probability weighted by the Boltzmann factor. We know though that, for example, when considering an asymmetric or a symmetric Eckart barrier and using the uniform semiclassical expression, the result is too low at high temperature and, in general, is within 10% of the exact answer.71 The steepest descent estimate is then within 20% of the exact answer and 10% of the numerically integrated uniform semiclassical answer. Moreover, the uniform semiclassical estimate does not give the correct leading order correction of order ℏ2. At high temperature (small u), it goes as159 
κuscβ1+u224118βV4V2253V32V23,
(4.23)
and this differs from the exact result given in Eq. (3.8).

In addition, employing the uniform theory at high temperatures is predicated on having a “good” expression for the action when the instanton energy increases above the barrier energy. This is not trivial if one wants to carry out the computations on the fly due to the difficulty in locating complex turning points. In this subsection, we describe some recent improvements, which answer these challenges.

An obvious failure of the uniform semiclassical expression for the energy dependent transmission coefficient [Eq. (2.10)] is that it always gives a probability of 1/2 when the action vanishes, that is, when the energy E = V. For example, for the symmetric Eckart barrier, at this energy, the exact transmission probability is
Tε=cosh8πα1cosh8πα+cosh8πα1α264,
(4.24)
where the dimensionless parameter α is proportional to ℏ. When α → 0, the probability is indeed ≃1/2, but as α grows, the transmission probability increases significantly at the barrier energy, as seen in Fig. 2. In typical computations, α ≃ 1. This is why in numerical studies the uniform semiclassical theory gives a rate estimate that is lower than the exact estimate.
FIG. 2.

The transmission probability at the barrier energy for a symmetric Eckart barrier is plotted as a function of the reduced ℏ parameter α=Md2ω. The exact (solid line), YF (dashed line), and V2 (dotted line) probabilities are shown. Note the good agreement especially when α is not much larger than unity.

FIG. 2.

The transmission probability at the barrier energy for a symmetric Eckart barrier is plotted as a function of the reduced ℏ parameter α=Md2ω. The exact (solid line), YF (dashed line), and V2 (dotted line) probabilities are shown. Note the good agreement especially when α is not much larger than unity.

Close modal
To correct for this, we note that without loss of generality, the exact energy dependent transmission probability for one-dimensional systems may always be written formally as in Eq. (2.10). The question is what is the exponent SE. Thus far, the only suggestion made was that it is the semiclassical Euclidean action, as derived from the semiclassical theory. However, there is nothing, in principle, that prevents us from using different expressions, provided that in the limit ℏ → 0, they reduce to the Euclidean action. For example, we may write the uniform transmission probability in two different ways,
TV2E=11+expSEΔE
(4.25)
or
TYFE=11+expSE+ΔS,
(4.26)
where the subscripts V2 and YF will be explained shortly. In the following, we will assume that the two constants ΔE and ΔS are independent of the energy and go as ℏ2. To determine them, we demand that in the high temperature limit, the resulting thermal rate gives the correct limit in terms of ℏ2 as given in Eq. (3.8). This also demonstrates the importance of deriving analytical results.
It is then a matter of straightforward expansion to find that159 
ΔE=2ω264V4V227V329V23
(4.27)
and
ΔS=2πΔEω.
(4.28)
Going back to the symmetric Eckart potential barrier, we note that if the fourth order derivative is positive, and indeed it is, then the energy and action shifts are negative, implying that the transmission probability will equal 1/2 at energies that are lower than the barrier height. The transmission probabilities at the barrier height for the symmetric Eckart barrier as functions of the reduced ℏ parameter [α=/Md2ω], where d is the length scale of the Eckart barrier obtained from the V2 and YF approximations as defined in Eqs. (4.25) and (4.26), respectively, using the results given in Eqs. (4.27) and (4.28), are compared with the exact results for the symmetric Eckart barrier in Fig. 2. One notes that foremost, the resulting shifts in the energy or action lead to an essentially exact result when α or equivalently ℏ is small. This is a reflection of the fact that the energy and action shifts were obtained by forcing the respective expressions to give the exact ℏ2 correction to the thermal transmission probability. Second, shifting the energy (V2) does not give as good agreement with the exact results as the constant action shift (YF), which is essentially exact up to rather large values of the reduced ℏ parameter α. This suggests that adding an ℏ2 correction term to the action appearing in the uniform semiclassical expression will give a more accurate approximation than obtained by shifting the energy.

It is due to the fact that the fourth derivative of the potential is positive, that the half point of the transmission probability is shifted to energies lower than the barrier energy, and thus, the thermal transmission coefficient is greater than unity. However, there are potentials for which the fourth derivative is negative and this could shift the midpoint to energies that are higher than the barrier energy and sometimes to thermal transmission probability ratios that are less than unity. In these cases, quantum above barrier reflection is the dominant quantum effect. Both the V2 and YF approaches can account for this quantum reflection as well. The CMD, RPMD, and classical Wigner approximations miss it.

We now relate to the notations V2 and YF. One of the central approximations used for obtaining thermal rates on the fly and at moderate to high temperatures is known as VPT2 theory.40,100 Its origins are in vibrational perturbation theory.38,39 Miller et al. used the well-known second order expansion of the energy in terms of the harmonic mode energies ωjnj+1/2 developed extensively in spectroscopy and inverted it for barrier crossing by considering one mode to have an imaginary action.32 This gave an analytic expression for the tunneling action in terms of the energy and up to the fourth derivatives of the potential at the barrier energy. Barker and co-workers then realized that one should also take into consideration that in the vibrational perturbation theory, there is a shift in the zero point energy, which is not accounted for by the classical Euclidean action.40,41 This shift is of order ℏ2 and affects the energy action relationship of VPT2 theory by shifting the energy with this same amount. It turns out that this energy shift in one dimension is precisely the energy shift parameter ΔE of Eq. (4.27). This is why VPT2 theory gives the correct ℏ2 term in the high temperature limit of the thermal rate, while CMD, RPMD, and classical Wigner do not. It is also the reason underlying the choice of V2 nomenclature used in Eq. (4.25).

There is though a subtle and meaningful difference between the approximation given in Eq. (4.25) and VPT2 theory. Due to the quadratic energy action relation, VPT2 theory forces a specific form for the dependence of the action on the energy. For the symmetric Eckart barrier, this relation is exact, but not so for other barriers, for example, the asymmetric Eckart barrier where VPT2 is known to fail in the deep tunneling region.44 If one uses the correct Euclidean action and only shifts its energy scale by ΔE, one significantly improves the VPT2 estimate for the rate also in the deep tunneling regime.159 The disadvantage is that, in general, one does not know the energy action relation, and this has to be established numerically.

There is a different way of getting good agreement between the uniform semiclassical theory and the exact results for asymmetric Eckart potentials. This was suggested by Eckart137 and implemented by Yasumori and Fueki,160 where each term in the exact transmission probability expression of the asymmetric Eckart barrier, which goes as coshx, is replaced by the positive exponential part expx/2. The resulting expression is almost identical to the uniform expression of Eq. (2.10) using the Euclidean action with one difference. The Euclidean action is shifted by a constant amount of order ℏ2. It is well known36,109,161 that for the asymmetric (and of course also the symmetric) Eckart barrier, as long as ℏ is not too large, the Yasumori–Fueki approximation is highly accurate. Its generalization to any one-dimensional system is given in Eq. (4.26), hence the notation YF. Can one do better? Perhaps yes, by merging the two forms, that is, by suggesting that
TE11+expSEΔE+ΔS,
(4.29)
retaining the action correction ΔS and the energy correction ΔE as parameters. Such a strategy has not yet been implemented as of the writing of this Perspective.

There remains a real challenge and that is how to generalize the one-dimensional results described in this subsection to multidimensional systems. This is the topic of Subsection IV D.

The semiclassical rate expression, based on vibrational perturbation theory derived by Miller et al.,32 is readily implemented on the fly. The central idea was to expand the Hamiltonian energy to second order in the harmonic energies for all degrees of freedom. Such expansions were considered to be “standard” in the vibrational spectroscopy community and led to rather accurate agreement between spectroscopic measurement and energy differences between vibrational energy levels, even for relatively highly excited vibrations.162 The assumption of Miller et al. was that if such a strategy works well for spectroscopy, it should do the same also for barrier crossing problems. The only difference would be that in barrier crossing, there is one negative normal mode eigenvalue at the saddle point so that the frequency and related coupling coefficients for this unstable mode would be imaginary. When considering N + 1 degrees of freedom, with the energy defined as zero at the bottom of the reactants well, the corresponding perturbation theory expression for the quantized energy, with normal mode frequencies ωk and anharmonic coupling coefficients xkk, was
En1,,nN+1=k=1N+1ωknk+12+kk=1N+1xkknk+12nk+12.
(4.30)
At a saddle point, one normal mode, denoted as F, has an imaginary frequency such that ωF=iωF.

The semiclassical quantization condition for a stable mode is well known to be SkE=2πnk+12 so that, in principle, one may identify each term in the Hamiltonian with a classical action. Miller et al.32 used this for only one mode—the unstable mode—with the identification nF+12=iθE/π, where 2θE=SE, the Euclidean action integral given in Eq. (2.11). Inserting this into the relation (4.30), one obtains a quadratic energy action relationship and the quadratic equation is readily solved. In this way, one obtains a tunneling action expression that depends on the quantum numbers of all the stable modes. Keeping the total energy E constant and summing over all states with total energy below E give a uniform semiclassical estimate for the cumulative transmission probability. The thermal rate is then obtained by appropriate Boltzmann averaging. The beauty of the theory is that all the coefficients in the energy action relation are expressed in terms of the second, third, and fourth order derivatives of the potential at the barrier.43 These derivatives are then computed from first-principles quantum chemistry programs, providing a relatively inexpensive ab initio based rate theory. In a way, it goes beyond instanton rate theory since it takes into consideration the nonlinear coupling between normal modes, while the instanton theory only has the linear terms related to the stability frequencies of the instanton.

As already mentioned, there was something missing. It was well understood within the spectroscopy community that the expansion appearing in Eq. (4.30) has an added constant term (ΔE), which to second order also depends on up to the fourth derivatives of the potential about the minimum. This constant contributes to the zero point energy but is not included in the expression given in Eq. (4.30) when all nk’s are set to 0. It was added to the VPT2-based rate theory by Barker and co-workers.40 We have already seen in Subsection IV C that this term is important as it correctly shifts the half point of the energy dependent transmission probability.

VPT2 theory is though incomplete, as we have seen. Especially for asymmetric systems and low temperatures, the dependence of the VPT2 action on the energy is too simplistic. It is also unreasonable to expect that one can derive a “good” approximation at all energies, including those much below the barrier energy, only based on a fourth order expansion of the potential about the barrier. A reasonable compromise would be to merge the VPT2 and instanton theories such that when the instanton energy is below the barrier, one uses the multidimensional microcanonical uniform instanton theory, while for energies above the barrier, one uses the VPT2 theory (a higher order approximation163 depends on up to the eighth order derivatives of the potential at the barrier, and as of the writing of this Perspective, these are too difficult to generate reliably from quantum chemistry codes). Such a methodology has many advantages. Coupling between modes becomes more important as the temperature is raised, and this could be well accounted for by VPT2 theory. At low temperature, the harmonic expansion implied by the instanton’s stability frequencies could in many cases be sufficient, especially for those modes whose ground state energy would be greater than kBT.

Finally, a comment on a multidimensional generalization of the YF approach. It can also be included in a VPT2-type framework. Instead of considering the energy relative to the multidimensional zero point energy contribution as in VPT2 theory and as shown in the one-dimensional V2 theory above [Eq. (4.25)], one could use the form of Miller et al. of the energy action relationship as in Eq. (4.30) without the energy shift, but include it in the definition of the action as in YF theory. Instead of the replacement nF+12=iSE/2π, one would use the YF relation nF+12=iSE+ΔS/2π, where the shift ΔS is related to the (multidimensional) zero point energy shift ΔE as in Eq. (4.28). This will not fully resolve the difficulties with VPT2 theory at low temperature, but may extend its validity to lower temperatures and more asymmetric systems. All this leads to the conclusion that there is work left for future investigators.

The difficulties associated with the “old” crossover temperature described above were not restricted to polyatomics. As already mentioned in the Introduction, the quantum thermal rate expression for dissipative systems derived by Wolynes,84 which is valid when the friction is moderate to strong—the spatial diffusion limited regime—and is the analog of the classical Kramers–Grote–Hynes expression,12,81,82 diverges at the “old” crossover temperature since it was derived on the basis of a parabolic barrier potential. It has been applied in the past to study experimental tunneling phenomena and isotope effects.164–166 It is therefore of practical interest to derive a similar expression, which does not have any divergence to it.

To reiterate briefly, dissipative systems are defined through a Hamiltonian in which the potential couples a “system” coordinate (q) bilinearly with a continuum of harmonic bath oscillators79 with mass weighted momenta and coordinates pxj,xj, frequencies ωj, and bilinear coupling coefficients cj,
H=pq22M+Vq+jpxj22+ωj22xjcjωj2Mq2.
(5.1)
When considered classically, it is well known that the classical equation of motion for the system coordinate may be recast as a generalized Langevin equation,80 in which the time dependent friction function is identified as
γt=jcj2ωj2cosωjt.
(5.2)
Kramers’ theory is aimed at assessing the rate of escape of a thermal particle, trapped in a well and escaping through a barrier. If the potential Vq is approximated as a parabolic barrier as in Eq. (3.1), the dissipative Hamiltonian has a quadratic form and so may be diagonalized using a normal mode transformation.83,108 If originally there are N harmonic bath modes, then this diagonalization gives N stable normal modes characterized by the frequencies λj and one unstable mode characterized by the barrier frequency λ, which turns out to be identical to the well-known Kramers–Grote–Hynes frequency whose relationship to the Laplace transform of the time dependent friction (denoted by a “hat”) is12,81,82
λ2+λγ̂λ=ω2.
(5.3)
Similarly, one may approximate the region of the well as a parabolic well so that again the dissipative Hamiltonian may be diagonalized, leading to N + 1 stable modes characterized by the frequencies λj. The escape rate is then obtained with the thermal equilibrium expression,85,
Γ=12πβZβZaβ=expβV2πβ×Πj=1Nsinhβλj2sinhβλ02sinhβλj2βλ2sinβλ2ω02πλωexpβVΞW,
(5.4)
where Zβ is the product of the partition functions of the N stable oscillators at the barrier with the “partition function” of the parabolic barrier, and similarly, Zaβ is the product of the partition functions of the N + 1 stable oscillators at the well. ω0 is the harmonic frequency at the well in the absence of friction. This result is identical to the Wolynes expression for the rate in the spatial diffusion limited regime.85 The last identity defines the Wolynes quantum factor,84 
ΞW=ωω0sinhβλ02sinβλ2πj=1Nsinhβλj2sinhβλj2,
(5.5)
which goes to unity in the limit that ℏ → 0.
Due to the parabolic barrier approximation, it is obvious that the Wolynes factor diverges at the friction dependent crossover temperature ℏβcλ = 2π, restricting its range of validity. One straightforward remedy would be to make use of the multi-dimensional instanton-based expression, as given in Eq. (4.20), identifying the stability frequencies of the instanton ωjE with the barrier stable normal mode frequencies λj. There are though a few difficulties. Foremost, one would have to find the instanton and its stability frequencies numerically. To do this, one needs to have a finite bath, and so the computation has to be converged with respect to the number of bath oscillators used. In contrast, the Wolynes factor may be expressed analytically in terms of the Laplace transform of the time dependent friction,85 
ΞW=πk=1ω02+νk2+νkγ̂νkω2+νk2+νkγ̂νk,
(5.6)
where νk=2πkβ is the k-th Matsubara frequency. This expression is readily converged with respect to the index k.
One notes though that the Wolynes factor may be rewritten as158,
ΞUniform=ΞWκpbβκpbβ=κpbβπk=1ω02+νk2+νkγ̂νkνk+λνk2νk+λ+γ̂νk,
(5.7)
where now the product of the Matsubara frequency terms no longer diverges, and the divergence is only in the parabolic barrier transmission factor. It is then very appealing to replace κpbβ with a non-parabolic barrier transmission factor, which does not diverge. This may be achieved in different ways, for example, a VPT2-based transmission factor. Or, one may use a separable approximation for the unstable mode motion and then obtain a one-dimensional instanton estimate of the thermal rate as above.158 The accuracy of such approaches remains to be studied by comparing with numerical benchmark quantum results. However, importantly, the Wolynes expression in its generalized form may be used without fear of having to deal with divergences.
The complete expression for the thermal rate of the turnover theory is87,88
Γβ=ω02πexpβVϒλωΞ,
(5.8)
where ϒ denotes the “depopulation factor,” which goes to unity when the friction becomes sufficiently strong. In Subsection V A, we discussed the spatial diffusion factor as expressed through the Wolynes ratio of partition functions Ξ. In this subsection, we discuss an instanton-based version of the depopulation factor.167 In the energy diffusion regime, the most important parameter is a reduced energy loss δ, which describes the energy lost by the particle as it traverses from the barrier to the repulsive wall of the well and back, or in the case of motion on a periodic potential, the loss of energy as the particle moves from one barrier to the adjacent one. When this energy loss is small, Kramers already showed that the rate goes as δ. When the energy loss is large, the depopulation factor goes to unity.

Originally, the Kramers turnover theory as developed by Mel’nikov and Meshkov86 and then formalized by PGH87 was classical. A short time later, Rips and Pollak88 developed a semiclassical version of the turnover theory, which was based on parabolic barrier transmission and reflection coefficients through the barrier. This theory by nature was constructed for temperatures above the parabolic barrier crossover temperature and so could not be extended to lower temperature. Here, we briefly describe how the uniform instanton theory is used to derive an instanton-based turnover theory, which does not suffer from divergences.

The central idea is to follow the previous derivation presented in Ref. 88, but instead of parabolic barrier reflection and transmission coefficients, one may use the uniform semiclassical theory. At a given temperature, one determines the relevant (reduced) instanton energy ɛβ = βEβ and corresponding action Sεβ. This action is then expanded linearly about the instanton energy,
SεSεβτεββεεβ,
(5.9)
where τεβ is the period of the instanton, which is of course proportional to the derivative of the action with respect to the energy. This form is analogous to the parabolic barrier form for the action, and therefore, the whole machinery used to develop the theory based on the parabolic barrier potential may be used smoothly to generalize the theory.
In words, in the instanton-based theory, the instanton trajectory is no longer at the barrier top, but at the instanton energy, which is below the barrier. The real time trajectory, which previously went from the barrier top to the turning point, is now a periodic trajectory going from the barrier region to the repulsive wall of the well and back. The energy loss is smaller, as it entails an orbit with less energy and a shorter period. The resulting expression is
ϒuniform=expβτεβsinπβτεβdy×ln1P̃yi2cosh2πβτεβycosπβτεβ,
(5.10)
where the double sided Laplace transform of the energy transfer Gaussian probability P̃yi2, which is central to the theory, is
P̃yi2=expδεβy2+14,
(5.11)
and it is here that the instanton energy dependent reduced energy loss δεβ appears.

In this way, an instanton-based semiclassical turnover theory has been formulated, which does away with the divergences resulting from parabolic barrier based theories. Is this the end? By no means. The semiclassical theory does not take into account the quantization of energy levels in the well. This implies that the last word has not been said at low temperature. This may be especially important when dealing with motion on a periodic potential, where on the one hand one will have bands of states whose width depends on the tunneling probability, and these bands in the absence of friction promote motion through many sites.

There is another aspect of the theory that one should remember; it is based on one-dimensional motion of the “system.” Its generalization to many degrees of freedom remains a challenge,99 although one could argue that in larger systems, the depopulation factor is generically unity due to the long time spent by the system in the multidimensional well, leading to extensive energy exchange as expressed through the multidimensional energy loss.168 

Finally, one notes that the instanton-based quantum version of turnover theory has not faced serious benchmark computations, which would check how accurate it is under various circumstances.

A central theme of this Perspective is to assess whether we have at our disposal today a reliable methodology for computing thermal reaction rates on the fly for chemical reactions, whether in the gas or condensed phases. It has been claimed62 that “Quantum mechanical effects on reaction rates can be enormous, even at room temperature, and the problem of including these effects in simulations of a wide variety of chemical reactions in complex systems has now effectively been solved.” Is this really the present status? In principle, there are numerically exact methods, such as the Multi-Configuration Time Dependent Hartree (MCTDH) method;46,47 however, it is limited, in practice, to systems with at most a few dozen (anharmonic) degrees of freedom. It is this limitation that forces us to consider approximate methods, some of which have been detailed above.

What do we mean when we say that the problem has been solved? Most approximate methods do not come with an objective criterion for the error made in the computation. Even if, in principle, one may create a perturbation series, which allows us to consider the leading order correction to the estimate, such as in semiclassical initial value representations151,169 or the classical Wigner methodology,49 in practice, such computations are prohibitive,33 especially on the fly. Artificial intelligence methods for obtaining “cost-effective” and reliable force fields are being developed170–173 so that maybe in the near future, one will find that the computation of correction terms will become feasible, but at present, such a program has not been implemented. Due to this limitation, we must consider other, less reliable, but “reasonable” criteria for reliability of the suggested solutions. Unfortunately, this injects a measure of subjectivity in the assessment. What may be considered by some as “reasonable” may be considered by others either as insufficient or conversely, demanding too much. Nevertheless, it is useful to discuss what may be expected so that in the least, when employing one of the methods, we have some understanding of its limitations and numerical accuracy.

Let us first consider the “high temperature limit.” In one dimension, in the limit u = ℏβω → 0, any method should reduce to the correct classical limit. However, what happens when considering more than one dimension? u may be sufficiently small such that the quantum tunneling or above barrier reflection correction may be ignored, yet at the saddle point, there may be modes whose frequency ωj is large in the sense that βωj1 and their motion is quantized. This reminds us that even in the high temperature limit, the method must account for quantization in the various degrees of freedom and this may be much more important than accounting for classical recrossing of a transition state. The CMD, RPMD, and classical Wigner dynamics methods do this correctly since, in principle, they use a numerically exact algorithm for computing the quantum thermal density. The VPT2 method accounts for this reasonably at high temperature, though to the best of my knowledge, a thorough comparison with the other three methods has not been implemented apart from the hydrogen exchange reaction where good agreement with numerically accurate computations has been demonstrated.174 We do not know how accurate the perturbation theory really is for molecular systems with many atoms, and this should be especially worrisome at high temperature, where inter-mode coupling becomes increasingly important. The “old” instanton method is not relevant in this context. The uniform instanton method described in this Perspective is relevant since it reduces at high temperature to the VPT2 form.

There is though another side to the high temperature limit. In one dimension, only the VPT2 theory and its YF variant are exact, giving the correct leading order term in ℏ2. The other theories are not. Is this a cause for real concern? Quantitatively probably not, as the differences between the quantum and classical rates in the high temperature parameter regimes are typically of the order of 20% or even less, and there are few experiments that will give rates with higher accuracy. Yet, this does demonstrate that there is something fundamental missing in some of the methods. Does this carry over to lower temperature regions? It is worrisome when one finds that some methods cannot distinguish whether quantum tunneling is more or less important than quantum reflection. Can one improve the approximations?

A real challenge is found when lowering the temperature to the region where u ≃ 2π. On the one hand, the “old” instanton does not exist in this region, or at best, when u is somewhat larger than 2π, it will give results that are too large due to the small denominator associated with the second derivative of the action, which is small in this temperature region. The formal justification for the RPMD or CMD methods in the sense that they reduce to the “old” instanton method at low temperature is not valid in this intermediate temperature regime. Can one further justify them by showing that they reduce to the uniform instanton result in this parameter region? This too is an open question. On the other hand, although the uniform instanton method does not suffer from the divergence of the “old” instanton method, all perpendicular degrees of freedom are treated linearly through the stability frequencies associated with the instanton. Is this accurate enough? We do not know, it is necessary to compare the multidimensional uniform theory with exact quantum benchmark results. In this context, one should also mention the development of machine learning algorithms for rate constants.175–177 One finds that precisely the “intermediate” temperature regime is the most challenging.177 

At low temperatures, in the regime of deep tunneling, there also remain open questions. In one dimension, we have seen that adding a term which goes as ℏ2 to the action improves its agreement with the exact quantum results significantly. It accounts correctly for the shift of the half point in the energy dependent transmission probability. The YF correction has thus far been tested only for one-dimensional systems. How does one generalize it? The resulting steepest descent rate estimate would be affected by this term, which is absent from the RPMD or CMD methods. One missing link comes from the high temperature regime, where one should derive analytically the leading order ℏ2 correction term to the rate in the multidimensional case and then be guided by the exact result to generalize the YF result.

One might argue that the future lies in machine learning based methods. Machine learning for force fields is a very active field.170–173 In the past few years, Bowman and co-workers have been using machine learning methods to predict reaction rates for systems with up to ten atoms.175–177 However, the quality of any machine learning method depends on its input. As long as the input is not precise, the resulting “predictions” will not be precise either. It is especially in this context that analytic results which are precise, are critical for assuring the quality of machine learning algorithms.

All these questions notwithstanding, it would seem that as of the writing of this Perspective, a “reasonable” middle way is to use the microcanonical uniform instanton YF theory in the range uπ and its VPT2 extension with the YF correction for uπ. From a numerical point of view, this is probably the cheapest possible theory. There is no need to average over trajectories, and such a methodology is readily amenable to on-the-fly computation. The energy averaging is not much of a restriction since typically the action and stability frequencies of the instanton change smoothly with energy. But is it sufficient? Not necessarily. At low temperature, none of these theories can account for either resonance phenomena or the quantum threshold reflection phenomenon.

One should add another proviso. The theories discussed thus far consider thermal or microcanonical transmission factors, with summation over all initial states. They do not give too much information about state to state reaction probabilities. Such information may be obtained with semiclassical initial value representation methods, but as already mentioned, these are as of now too expensive for on-the-fly implementation. This whole discussion can be thus boiled down to one sentence—quantum rate theory will continue to challenge us at least in the near future.

The discussion in this Perspective was limited to motion on a single Born–Oppenheimer potential energy surface. Using the imaginary free energy method,123 it was extended to non-adiabatic transitions by Cao and Voth.178,179 The same challenge has been treated within a Fermi golden rule framework in Ref. 180. Can one use the uniform instanton theory also in this context? Probably yes, but this needs to be implemented. There are other important issues missing. Reactions in liquids present serious challenges, which at this point are even difficult to handle on the fly within the realm of classical mechanics.181 

A second aspect of this Perspective was a discussion of rate theory in dissipative systems. Why was this included here? For a number of reasons. As we have already seen, analytic theories help us understand and analyze the relevant quantum phenomena. Removal of the divergence in the spatial diffusion limited regime of dissipative systems implies that the theory is more amenable for analysis of experimental results. For example, tunneling has been established for H and D atom diffusion on metal surfaces,164 but a good theoretical interpretation is lacking. How do phonons affect the rate? Is electronic friction important?182 Progress made in understanding the uniform instanton theory has led to an analytic instanton-based turnover theory; it may now be used to answer the following questions: How does friction affect the hopping distribution over different wells? How important is above barrier reflection when determining diffusion coefficients? The theory does have limitations. It is presented in terms of memory friction, but the memory time must not be longer than the period of motion in the well at the instanton energy.183 What happens at longer memory times? At low energies, especially when the friction is weak, quantization of energy levels in the well of the reacting system should play a role, and this is not included in the theory. The present turnover theory has not been challenged with numerically exact computations. The RPMD method has been used to study the quantum turnover problem.182,184 Does the present analytic turnover theory agree with the RPMD results? The quantum turnover theory was originally devised for the study of Josephson junction circuits.4,79,89 Does the present theory shed any further light on these systems or will it lead to only marginally corrected results?

Hopefully, this rather personal perspective will establish once again that although numerics are important, there is place for qualitative ideas expressed through analytic theories even in this day and age. Tolman’s remarks of almost one hundred years ago on reaction velocities are still valid today.

I thank Professor Michele Ceotto for his detailed comments on an earlier version of this perspective, which significantly improved it. I also thank Professor Jianshu Cao and Professor Jian Liu and my Ph.D. student Mr. Sameernandan Upadhyayula for useful discussions. This work was graciously supported by a grant from the Israel Science Foundation.

The author has no conflicts to disclose.

Eli Pollak: Conceptualization (equal); Formal analysis (equal); Funding acquisition (equal); Investigation (equal); Methodology (equal); Validation (equal); Writing – original draft (equal); Writing – review & editing (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.

1.
2.
K. J.
Laidler
,
J. Chem. Educ.
61
,
494
(
1984
).
3.
R. C.
Tolman
,
Statistical Mechanics with Applications to Physics and Chemistry
,
ACS Monograph Series
(
The Chemical Catalog Co.
,
New York
,
1927
), p.
323
.
4.
P.
Hänggi
,
P.
Talkner
, and
M.
Borkovec
,
Rev. Mod. Phys.
62
,
251
(
1990
).
5.
6.
H.
Pelzer
and
E. P.
Wigner
,
Z. Phys. Chem.
15B
,
445
(
1932
).
7.
E. P.
Wigner
,
Z. Phys. Chem.
19B
,
203
(
1932
).
8.
E. P.
Wigner
,
Trans. Faraday Soc.
34
,
29
(
1938
).
9.
H.
Eyring
,
J. Chem. Phys.
3
,
107
(
1935
).
10.
H.
Eyring
,
Trans. Faraday Soc.
34
,
41
(
1938
).
11.
E.
Pollak
and
P.
Talkner
,
Chaos
15
,
026116
(
2005
).
13.
The Transition State: A Symposium Held at Sheffield on 3rd and 4th April 1962, Chemical Society Great Britain, 1962. Discussion on the Paper by H. Eyring et al., p.
25
.
14.
D.
Chandler
,
J. Chem. Phys.
68
,
2959
(
1978
).
15.
R. P.
Bell
,
The Tunnel Effect in Chemistry
(
Chapman & Hall
,
London
,
1980
).
16.
Tunneling
, edited by
J.
Jortner
and
B.
Pullman
(
D. Reidel Publ. Co.
,
1986
).
17.
V. I.
Gol’danskii
,
L. I.
Trakhtenberg
, and
V. N.
Fleurov
,
Tunneling Phenomena in Chemical Systems
(
Gordon & Breach
,
1988
).
18.
Chemical Dynamics at Low Temperatures
,
Series Advances in Chemical Physics Vol. 88
, edited by
V. A.
Benderskii
,
D. E.
Makarov
, and
C. A.
Wight
(
Wiley
,
New York
,
1994
) and references therein.
19.
Tunneling and its Implications
, edited by
D.
Mugnai
,
A.
Ranfagni
, and
L. S.
Schulman
(
World Scientific
,
1997
).
20.
J.
Ankerhold
, in
Quantum Tunneling in Complex Systems: The Semiclassical Approach
,
Springer Tracts in Modern Physics Vol. 224
(
Springer Verlag
, Berlin,
2007
).
21.
M.
Razavy
,
Quantum Theory of Tunneling
, 2nd ed. (
World Scientific
,
2013
).
22.
Tunnelling in Molecules
, edited by
J.
Kästner
and
S.
Kozuch
(
Royal Society of Chemistry
,
2012
).
23.
R. D.
Levine
and
S. F.
Wu
,
Chem. Phys. Lett.
11
,
557
(
1971
).
24.
J. M.
Bowman
and
A.
Kuppermann
,
Chem. Phys. Lett.
12
,
1
(
1971
).
25.
A. B.
Henson
,
S.
Gersten
,
Y.
Shagam
,
J.
Narevicius
, and
E.
Narevicius
,
Science
338
,
234
(
2012
).
26.
P.
Pechukas
, in
Dynamics of Molecular Collision B
, edited by
W. H.
Miller
(
Springer
,
1976
), Chap. 6.
27.
W. H.
Miller
,
J. Chem. Phys.
61
,
1823
(
1974
).
28.
F. J.
McLafferty
and
P.
Pechukas
,
Chem. Phys. Lett.
27
,
511
(
1974
).
29.
W. H.
Miller
,
S. D.
Schwartz
, and
J.
Tromp
,
J. Chem. Phys.
79
,
4889
(
1983
).
30.
M.
Topaler
and
N.
Makri
,
J. Chem. Phys.
101
,
7500
(
1994
).
31.
W. H.
Miller
,
J. Chem. Phys.
62
,
1899
(
1975
).
32.
W. H.
Miller
,
R.
Hernandez
,
N. C.
Handy
,
D.
Jayatilaka
, and
A.
Willetts
,
Chem. Phys. Lett.
172
,
62
(
1990
).
33.
E.
Martin-Fierro
and
E.
Pollak
,
J. Chem. Phys.
126
,
164108
(
2007
).
34.
R. P.
Bell
,
Trans. Faraday Soc.
55
,
1
(
1959
).
35.
E.
Pollak
and
J.
Cao
,
J. Chem. Phys.
157
,
074109
(
2022
).
36.
E.
Pollak
and
J.
Cao
,
Phys. Rev. A
107
,
022203
(
2023
);
E.
Pollak
and
J.
Cao
,
Phys. Rev. A
109
,
039901
(
2024
).
37.
R.
Hernandez
and
W. H.
Miller
,
Chem. Phys. Lett.
214
,
129
(
1993
).
38.
39.
S.
Califano
,
Vibrational States
(
John Wiley & Sons
,
1976
), Chap. 9.
40.
T. L.
Nguyen
,
J. F.
Stanton
, and
J. R.
Barker
,
Chem. Phys. Lett.
499
,
9
(
2010
).
41.
T. L.
Nguyen
,
J. F.
Stanton
, and
J. R.
Barker
,
J. Phys. Chem. A
115
,
5118
(
2011
).
42.
C.
Aieta
,
F.
Gabas
, and
M.
Ceotto
,
J. Chem. Theory Comput.
15
,
2142
(
2019
).
43.
G.
Mandelli
,
L.
Corneo
, and
C.
Aieta
,
J. Phys. Chem. Lett.
14
,
9996
(
2023
).
44.
P.
Goel
and
J. F.
Stanton
,
J. Chem. Phys.
149
,
134109
(
2018
).
45.
A. F.
Wagner
,
J. Phys. Chem. A
117
,
13089
(
2013
).
46.
H.-D.
Meyer
,
U.
Manthe
, and
L.
Cederbaum
,
Chem. Phys. Lett.
165
,
73
(
1990
).
47.
M. H.
Beck
,
A.
Jäckle
,
G. A.
Worth
, and
H.-D.
Meyer
,
Phys. Rep.
324
,
1
(
2000
).
48.
E. J.
Heller
,
J. Chem. Phys.
62
,
1544
(
1975
).
49.
J.
Shao
,
J.-L.
Liao
, and
E.
Pollak
,
J. Chem. Phys.
108
,
9711
(
1998
).
50.
J.
Liu
and
W. H.
Miller
,
J. Chem. Phys.
131
,
074113
(
2009
).
51.
J.
Liu
,
Int. J. Quantum Chem.
115
,
657
670
(
2015
).
52.
W. H.
Miller
,
J. Phys. Chem. A
105
,
2942
(
2001
).
53.
J.
Tatchen
and
E.
Pollak
,
J. Chem. Phys.
130
,
041103
(
2009
).
54.
M.
Ceotto
,
S.
Atahan
,
S.
Shim
,
G.
Tantardini
, and
A.
Aspuru-Guzik
,
Phys. Chem. Chem. Phys.
11
,
3861
(
2009
).
55.
C.
Mollica
and
J.
Vanicek
,
Phys. Rev. Lett.
107
,
214101
(
2011
).
56.
M.
Wehrle
,
M.
Šulc
, and
J.
Vaníek
,
J. Chem. Phys.
140
,
244114
(
2014
).
57.
M.
Buchholz
,
E.
Fallacara
,
F.
Gottwald
,
M.
Ceotto
,
F.
Grossmann
, and
S. D.
Ivanov
,
Chem. Phys.
515
,
231
(
2018
).
58.
J.
Cao
and
G. A.
Voth
,
J. Chem. Phys.
100
,
5106
5117
(
1994
).
59.
G. A.
Voth
, in
Progress in Theoretical Chemistry and Physics
, edited by
S. D.
Schwartz
(
Kluwer
,
Dordrecht
,
2000
), Vol.
5
, pp.
47
68
.
60.
S.
Habershon
,
D. E.
Manolopoulos
,
T. E.
Markland
, and
T. F.
Miller
III
,
Annu. Rev. Phys. Chem.
64
,
387
413
(
2013
).
61.
Y. V.
Suleimanov
,
F. J.
Aoiz
, and
H.
Guo
,
J. Phys. Chem. A
120
(
43
),
8488
8502
(
2016
).
62.
J. E.
Lawrence
and
D. E.
Manolopoulos
,
Faraday Discuss.
221
,
9
(
2020
).
63.
S.
Coleman
,
Phys. Rev. D
15
,
2929
2936
(
1977
).
64.
M.
Kryvohuz
,
J. Chem. Phys.
134
,
114103
(
2011
).
65.
M.
Kryvohuz
and
R. A.
Marcus
,
J. Chem. Phys.
137
,
234304
(
2012
).
67.
J.
Meisner
and
J.
Kästner
,
Angew. Chem., Int. Ed.
55
,
5400
5413
(
2016
).
68.
A. N.
Beyer
,
J. O.
Richardson
,
P. J.
Knowles
,
J.
Rommel
, and
S. C.
Althorpe
,
J. Phys. Chem. Lett.
7
,
4374
4379
(
2016
).
69.
J. O.
Richardson
,
Faraday Discuss.
195
,
49
(
2016
).
70.
J. O.
Richardson
,
Int. Rev. Phys. Chem.
37
,
171
(
2018
).
71.
S.
Upadhyayula
and
E.
Pollak
,
J. Phys. Chem. Lett.
14
,
9892
(
2023
).
72.
D. V.
Shalashilin
and
M. S.
Child
,
J. Chem. Phys.
115
,
5367
(
2001
).
73.
W. H.
Miller
,
Y.
Zhao
,
M.
Ceotto
, and
S.
Yang
,
J. Chem. Phys.
119
,
1329
(
2003
).
74.
M.
Ceotto
and
W. H.
Miller
,
J. Chem. Phys.
120
,
6356
(
2004
).
75.
J.
Vanicek
,
Rate Constant Calculation for Thermal Reactions: Methods and Applications
(
Wiley
,
Hoboken, NJ
,
2012
), pp.
67
92
.
76.
L.
Vaillant
,
M. J.
Thapa
,
J.
Vaníček
, and
J. O.
Richardson
,
J. Chem. Phys.
151
,
144111
(
2019
).
77.
R.
Iftimie
,
P.
Minary
, and
M. E.
Tuckerman
,
Proc. Natl. Acad. Sci. U. S. A.
102
,
6654
(
2005
).
78.
D.
Marx
and
J.
Hutter
,
Ab Initio Molecular Dynamics: Basic Theory and Advanced Methods
(
Cambridge University Press
,
2010
).
79.
U.
Weiss
,
Quantum Dissipative Systems
, 5th ed. (
World Scientific Publ. Co.
,
2021
).
80.
R.
Zwanzig
,
J. Stat. Phys.
9
,
215
(
1973
).
81.
R. F.
Grote
and
J. T.
Hynes
,
J. Chem. Phys.
73
,
2715
(
1980
).
82.
P.
Hänggi
and
F.
Mojtabai
,
Phys. Rev. A
26
,
1168
(
1982
).
83.
E.
Pollak
,
J. Chem. Phys.
85
,
865
(
1986
).
84.
P. G.
Wolynes
,
Phys. Rev. Lett.
47
,
968
(
1981
).
86.
V. I.
Mel’nikov
and
S. V.
Meshkov
,
J. Chem. Phys.
85
,
1018
(
1986
).
87.
E.
Pollak
,
H.
Grabert
, and
P.
Hänggi
,
J. Chem. Phys.
91
,
4073
(
1989
).
88.
I.
Rips
and
E.
Pollak
,
Phys. Rev. A
41
,
5366
(
1990
).
90.
L.
Rondin
,
J.
Gieseler
,
F.
Ricci
,
R.
Quidant
,
C.
Dellago
, and
L.
Novotny
,
Nat. Nanotechnol.
12
,
1130
(
2017
).
91.
P. L. G.
Müller
,
R.
Hernandez
,
R. M.
Benito
, and
F.
Borondo
,
J. Chem. Phys.
137
,
204301
(
2012
).
92.
Y.
Shigemitsu
and
Y.
Ohga
,
J. Solution Chem.
43
,
1746
(
2014
).
93.
B. R.
Ferrer
,
J. R.
Gomez-Solano
, and
A. V.
Arzola
,
Phys. Rev. Lett.
126
,
108001
(
2021
).
94.
A.
Militaru
,
M.
Innerbichler
,
M.
Frimmer
,
F.
Tebben- johanns
,
L.
Novotny
, and
C.
Dellago
,
Nat. Commun.
12
,
2446
(
2021
).
95.
F.
Ginot
,
J.
Caspers
,
M.
Krüger
, and
C.
Bechinger
,
Phys. Rev. Lett.
128
,
028001
(
2022
).
96.
D. C.
Senft
and
G.
Ehrlich
,
Phys. Rev. Lett.
74
,
294
(
1995
).
97.
T. R.
Linderoth
,
S.
Horch
,
E.
Lægsgaard
,
I.
Stens- gaard
, and
F.
Besenbacher
,
Phys. Rev. Lett.
78
,
4978
(
1997
).
98.
J.
Jacobsen
,
K. W.
Jacobsen
, and
J. P.
Sethna
,
Phys. Rev. Lett.
79
,
2843
(
1997
).
99.
E.
Hershkovitz
,
P.
Talkner
,
E.
Pollak
, and
Y.
Georgievskii
,
Surf. Sci.
421
,
73
(
1999
).
100.
X.
Shan
,
T. A. H.
Burd
, and
D. C.
Clary
,
J. Phys. Chem. A
123
,
4639
(
2019
).
101.
E.
Pollak
and
R.
Ianconescu
,
J. Phys. Chem. A
120
,
3155
(
2016
).
102.
R.
Ianconescu
and
E.
Pollak
,
J. Chem. Phys.
151
,
024703
(
2019
).
103.
E.
Pollak
and
S.
Miret-Artés
,
ChemPhysChem
24
,
e202300272
(
2023
).
104.
J.
Costley
and
P.
Pechuka
,
Chem. Phys. Lett.
83
,
139
(
1981
).
105.
R.
Kubo
,
M.
Yokota
, and
S.
Nakajima
,
J. Phys. Soc. Jpn.
12
,
1203
(
1957
).
106.
G. A.
Voth
,
D.
Chandler
, and
W. H.
Miller
,
J. Phys. Chem.
93
,
7009
(
1989
).
107.
X.
Sun
,
H.
Wang
, and
W. H.
Miller
,
J. Chem. Phys.
109
,
4190
(
1998
).
108.
J. L.
Liao
and
E.
Pollak
,
Chem. Phys.
268
,
295
(
2001
).
109.
E.
Pollak
,
S.
Upadhyayula
, and
J.
Liu
,
J. Chem. Phys.
156
,
244101
(
2022
).
110.
E. P.
Wigner
,
J. Chem. Phys.
5
,
720
(
1937
).
111.
112.
T.
Komatsuzaki
and
R. S.
Berry
,
J. Mol. Struct.: THEOCHEM
506
,
55
(
2000
).
113.
S.
Wiggins
,
L.
Wiesenfeld
,
C.
Jaffé
, and
T.
Uzer
,
Phys. Rev. Lett.
86
,
5478
(
2001
).
114.
T.
Bartsch
,
R.
Hernandez
, and
T.
Uzer
,
Phys. Rev. Lett.
95
,
058301
(
2005
).
115.
Y.
Nagahata
,
R.
Hernandez
, and
T.
Komatsuzaki
,
J. Chem. Phys.
155
,
210901
(
2021
).
116.
E.
Pollak
,
J. Chem. Phys.
74
,
6765
(
1981
).
117.
E.
Pollak
and
D.
Proselkov
,
Chem. Phys.
170
,
265
(
1993
).
119.
M.
Hillery
,
R. F.
O’Connell
,
M. O.
Scully
, and
E. P.
Wigner
,
Phys. Rep.
106
,
121
(
1984
).
121.
N.
Fröman
and
P. O.
Fröman
,
JWKB Approximation: Contributions to the Theory
(
North-Holland
,
Amsterdam
,
1965
).
122.
C. G.
Callan
and
S.
Coleman
,
Phys. Rev. D
16
,
1762
(
1977
).
124.
J. O.
Richardson
,
J. Chem. Phys.
144
,
114106
(
2016
).
125.
J. E.
Lawrence
and
J. O.
Richardson
,
Faraday Discuss.
238
,
204
(
2022
).
126.
J.
Cao
and
G.
Voth
,
J. Chem. Phys.
105
,
6856
(
1996
).
127.
128.
J.-L.
Liao
and
E.
Pollak
,
J. Chem. Phys.
111
,
7244
(
1999
).
129.
Y.
Zheng
and
E.
Pollak
,
J. Chem. Phys.
114
,
9741
(
2001
).
130.
J.-L.
Liao
and
E.
Pollak
,
J. Chem. Phys.
116
,
2718
(
2002
).
131.
Y.
Zheng
,
J. Chem. Phys.
122
,
094316
(
2005
).
132.
J.
Ankerhold
,
M.
Saltzer
, and
E.
Pollak
,
J. Chem. Phys.
116
,
5925
(
2002
).
133.
I. R.
Craig
and
D. E.
Manolopoulos
,
J. Chem. Phys.
122
,
084106
(
2005
).
134.
I. R.
Craig
and
D. E.
Manolopoulos
,
J. Chem. Phys.
123
,
034102
(
2005
).
135.
C.
Doubleday
,
R.
Armas
,
D.
Walker
,
C. V.
Cosgriff
, and
E. M.
Greer
,
Angew. Chem.
129
,
13279
(
2017
).
136.
C.
Castro
and
W. L.
Karney
,
Angew. Chem., Int. Ed.
59
,
8355
(
2020
).
140.
S.
Jang
,
A. V.
Sinitskiy
, and
G. A.
Voth
,
J. Chem. Phys.
140
,
154103
(
2014
).
141.
S.
Jang
and
G. A.
Voth
,
J. Chem. Phys.
144
,
084110
(
2016
).
142.
S.
Jang
and
G. A.
Voth
,
J. Chem. Phys.
146
,
174106
(
2017
).
143.
T. J. H.
Hele
and
S. C.
Althorpe
,
J. Chem. Phys.
138
,
084108
(
2013
).
144.
T. J. H.
Hele
and
S. C.
Althorpe
,
J. Chem. Phys.
139
,
084115
(
2013
).
145.
T. J. H.
Hele
and
S. C.
Althorpe
,
J. Chem. Phys.
144
,
174107
(
2016
).
146.
147.
P.
Hänggi
and
W.
Hontscha
,
J. Chem. Phys.
88
,
4094
(
1988
).
148.
J. B.
Rommel
,
T. P. M.
Goumans
, and
J.
Kästner
,
J. Chem. Theory Comput.
7
,
690
698
(
2011
).
149.
T.
Lamberts
,
P. K.
Samanta
,
A.
Köhn
, and
J.
Kästner
,
Phys. Chem. Chem. Phys.
18
,
33021
(
2016
).
150.
E.
Han
,
W.
Fang
,
M.
Stamatakis
,
J. O.
Richardson
, and
J.
Chen
,
J. Phys. Chem. Lett.
13
,
3173
3181
(
2022
).
151.
E.
Pollak
and
J.
Shao
,
J. Phys. Chem. A
107
,
7112
(
2003
).
152.
S.
Zhang
and
E.
Pollak
,
Phys. Rev. Lett.
91
,
190201
(
2003
).
153.
H.
Cartarius
and
E.
Pollak
,
Chem. Phys.
399
,
135
(
2012
).
154.
J. E.
Lawrence
,
J.
Dusek
, and
J. O.
Richardson
,
J. Chem. Phys.
159
,
014111
(
2023
).
155.
S.
Mandrà
,
S.
Valleau
, and
M.
Ceotto
,
Int. J. Quantum Chem.
113
,
1722
(
2013
).
156.
J. E.
Lennard-Jones
and
A. F.
Devonshire
,
156
,
6
(
1936
).
157.
158.
E.
Pollak
,
J. Chem. Phys.
159
,
224107
(
2023
).
159.
S.
Upadhyayula
and
E.
Pollak
,
J. Phys. Chem. A
(submitted).
160.
I.
Yasumori
and
K.
Fueki
,
J. Chem. Phys.
22
,
1790
(
1954
).
161.
I.
Shavitt
,
J. Chem. Phys.
31
,
1359
(
1959
).
162.
Q.
Yu
and
J. M.
Bowman
,
Mol. Phys.
113
,
3964
(
2015
).
163.
J.
Stanton
,
J. Phys. Chem. Lett.
7
,
2708
(
2016
).
164.
A. P.
Jardine
,
E. Y. M.
Lee
,
D. J.
Ward
,
G.
Alexandrowicz
,
H.
Hedgeland
,
W.
Allison
,
J.
Ellis
, and
E.
Pollak
,
Phys. Rev. Lett.
105
,
136101
(
2010
).
165.
E.
Pollak
,
J. Phys. Chem. B
116
,
12966
(
2012
).
166.
Y.
Liu
,
Y.
Yan
,
T.
Xing
, and
Q.
Shi
,
J. Phys. Chem. B
125
,
5959
(
2021
).
168.
E.
Hershkovitz
and
E.
Pollak
,
J. Chem. Phys.
106
,
7678
(
1997
).
169.
E.
Pollak
,
Springer Ser. Chem. Phys.
83
,
259
(
2007
).
170.
C.
Schran
,
J.
Behler
, and
D.
Marx
,
J. Chem. Theory Comput.
16
,
88
(
2020
).
171.
C.
Qu
,
P. L.
Houston
,
R.
Conte
,
A.
Nandi
, and
J.
Bowman
,
J. Phys. Chem. Lett.
12
,
4902
(
2021
).
172.
M. S.
Chen
,
J.
Lee
,
H.-Z.
Ye
,
T. C.
Berkelbach
,
D. R.
Reichman
, and
T. E.
Markland
,
J. Chem. Theory Comput.
19
,
4510
(
2023
).
173.
K.
Brezina
,
H.
Beck
, and
O.
Marsalek
,
J. Chem. Theory Comput.
19
,
6589
(
2023
).
174.
T. L.
Nguyen
,
J. R.
Barker
, and
J. F.
Stanton
,
Adv. Atmos. Chem.
1
,
403
492
(
2017
).
175.
P. L.
Houston
,
A.
Nandi
, and
J. M.
Bowman
,
J. Phys. Chem. Lett.
10
,
5250
(
2019
).
176.
A.
Nandi
,
J. M.
Bowman
, and
P. A.
Houston
,
J. Phys. Chem. A
124
,
5746
(
2020
).
177.
P. L.
Houston
,
A.
Nandi
, and
J. M.
Bowman
,
J. Phys. Chem. A
126
,
5672
(
2022
).
178.
J.
Cao
and
G. A.
Voth
,
J. Chem. Phys.
106
,
1769
(
1997
).
179.
S.
Jang
and
J.
Cao
,
J. Chem. Phys.
114
,
9959
(
2001
).
180.
M. J.
Thapa
,
W.
Fang
, and
J. O.
Richardson
,
J. Chem. Phys.
150
,
104107
(
2019
).
181.
E.
Grifoni
,
G.
Piccini
, and
M.
Parrinello
,
Proc. Natl. Acad. Sci. U. S. A.
116
,
4054
(
2019
).
182.
Y.
Litman
,
E. S.
Pós
,
C. L.
Box
,
R.
Martinazzo
,
R. J.
Maurer
, and
M.
Rossi
,
J. Chem. Phys.
156
,
194107
(
2022
).
183.
R.
Ianconescu
and
E.
Pollak
,
J. Chem. Phys.
143
,
104104
(
2015
).
184.
Y.
Litman
,
E. S.
Pós
,
C. L.
Box
,
R.
Martinazzo
,
R. J.
Maurer
, and
M.
Rossi
,
J. Chem. Phys.
156
,
194106
(
2022
).