We explore the relation between the quantum and semiclassical instanton approximations for the reaction rate constant. From the quantum instanton expression, we analyze the contributions to the rate constant in terms of minimum-action paths and find that two such paths dominate the expression. For symmetric barriers, these two paths join together to describe the semiclassical instanton periodic orbit. However, for asymmetric barriers, one of the two paths takes an unphysically low energy and dominates the expression, leading to order-of-magnitude errors in the rate predictions. Nevertheless, semiclassical instanton theory remains accurate. We conclude that semiclassical instanton theory can be obtained directly from the semiclassical limit of the quantum instanton only for symmetric systems. We suggest a modification of the quantum instanton approach which avoids sampling the spurious path and thus has a stronger connection to semiclassical instanton theory, giving numerically accurate predictions even for very asymmetric systems in the low temperature limit.

## I. INTRODUCTION

The inclusion of nuclear quantum effects in molecular dynamics calculations is a challenging task. One of the main difficulties in calculating exact quantum dynamics is the need to solve the time-dependent Schrödinger equation with many degrees of freedom. For atomistic descriptions of most molecular systems, the dynamics are (currently) computationally impossible to solve exactly due to the exponentially increasing size of the basis sets required for these calculations. For this reason, approximate quantum dynamics methods based on path integrals are increasingly favored.

A rigorous formulation of reaction rate theory is provided by the flux-side or flux-flux time correlation function formalism.^{1} These correlation functions can be evaluated with approximate quantum-dynamics methods, including linearized semiclassical initial-value representation (LSC-IVR),^{2} centroid molecular dynamics,^{3,4} ring-polymer molecular dynamics (RPMD),^{5–7} and the recently developed Matsubara dynamics^{8,9} and its approximations.^{10,11} The main approximation made by these methods is to neglect quantum coherences, which is not expected to be a cause of significant error for complex systems in thermodynamic equilibrium.^{2} Thus, solving the system exactly to obtain the long-time dynamics is, in many cases, unnecessary as the relevant processes can be approximated to a good accuracy at short times. In fact, if the problem is well defined, no time propagation is necessary at all and one can turn to so-called quantum transition state theories (QTSTs), where the expressions involving time correlation functions are approximated using only their properties at zero time.^{12,13}

Semiclassical instanton (SCI) theory provides one of the simplest formulations of a QTST for multidimensional tunneling.^{14,15} It can be derived rigorously by taking asymptotic (*ℏ* → 0) approximations to the flux-flux correlation function^{16} and is known in a number of formulations,^{17} including one obtained from the “ImF” approximation,^{18} which can all be shown to be equivalent.^{19} SCI theory is defined in terms of the dominant tunneling pathway located on the full-dimensional potential-energy surface and employs a harmonic approximation for fluctuations around this pathway. The ring-polymer instanton approach^{20–22} provides a computationally efficient algorithm for applying the theory in practice and can be automated for application to complex systems.^{23}

Although SCI theory is very powerful, the harmonic expansion about a single dominant tunneling path can introduce errors, especially for low-frequency anharmonic systems. An alternative method, known as the quantum instanton (QI),^{24} was developed as an attempt to correct the quantitative deficiencies of the semiclassical approximation. This method goes beyond the harmonic approximation employed by SCI theory and samples paths using efficient path-integral Monte Carlo or molecular dynamics methods.^{25} The QI method has been used to calculate reaction rates in many multidimensional systems^{26–28} and is a particularly efficient approach for the direct calculation of kinetic isotope effects.^{29–32} However, for large asymmetries (meaning the reaction is highly exothermic or endothermic), the rate constant predicted by the QI method has been observed to give large errors, whereas SCI theory remains well-behaved.^{17,33}

There are also a number of other path-integral-based QTSTs developed to go beyond the SCI approximation. Unlike in classical mechanics,^{34} it is not possible to derive them directly from the flux-side correlation function, as this function tends to zero at zero time if used in its standard form. Instead, QTSTs have been proposed based on a connection to semiclassical instanton theory, including a free-energy version of instanton theory^{35} and ring-polymer transition-state theory, which is itself related to the RPMD rate theory.^{21} More recently, ring-polymer transition-state theory has been rederived from a generalized flux-side correlation function, which yields a good approximation to the rate at zero time.^{36,37} This provides an extension of the centroid-based QTST of Voth, Chandler, and Miller,^{38} a method which does not dominantly sample the instanton and thus fails for asymmetric barriers.^{39}

In summary, it appears that the QTSTs which work best are those with the strongest connection to the SCI. In this article, we examine the relation between the QI and SCI theories in order to investigate the very plausible conjecture that the SCI is a semiclassical/steepest-descent approximation to the QI; put another way, the QI method is an improvement on the SCI with more accurate sampling of the paths. Surprisingly, we find that the two methods are not directly connected to each other except in the case of symmetric systems. We perform a semiclassical analysis of the quantities used in the QI method and find that the dominant paths that contribute are not the same as those which define the instanton periodic orbit. This leads to the observed error of the QI method for systems with large asymmetries and low temperatures. We then suggest an approach for modifying the QI expression to enforce sampling of the instanton orbit, which is seen to greatly improve the results.

This paper is organized as follows. We first give a derivation of the QI method from first principles in Sec. II, followed by a short investigation of the breakdown of the predicted rate for asymmetric systems. Motivated by this breakdown, we analyze the contributions from semiclassical paths in Sec. III, and we show that the SCI expression can only be recovered if the spurious paths are removed. We suggest adding a projection operator to solve this problem and describe an improved approach, which we call the projected quantum instanton method, in Sec. IV. We follow this with a discussion of the numerical results of the various methods in Sec. V, including the semiclassical approximations, which justifies the analysis in terms of semiclassical paths, before concluding in Sec. VI.

## II. QUANTUM INSTANTON APPROXIMATION

In Sec. II A, we re-examine the derivation of the QI method for calculating reaction rates. The derivation is based on a saddle-point approximation to the time integral of the flux-flux correlation function, and we show that the method breaks down for asymmetric barriers at low temperatures. The goal is to have a self-contained discussion of the QI method before undertaking a more detailed semiclassical analysis. There are two versions of the QI method: the original, which we shall continue to call QI, and a variant, which we call the second-order cumulant expansion (2OCE).

Although we limit ourselves to one dimension throughout the paper for simplicity, multidimensional extensions of the QI method exist.^{24–26,29} Our conclusions remain valid for these multidimensional cases.

### A. Derivation

We begin by rederiving the QI approximation from the Miller–Schwartz–Tromp expression for the reaction rate constant,^{1}

Here, *Q*_{r} is the reactant partition function and *C*_{ff} is the symmetrized flux-flux time correlation function,

generalized to two dividing surfaces. The Hamiltonian for a particle in a one-dimensional potential, *V* (*x*), is $\u0124=p^2/2m+V(x^)$, and the flux operators are given by

where *m* is the mass and $\delta ^j=\delta (x^\u2212xj)$ indicates a Dirac delta function centered at the *j*th dividing surface located at *x* = *x*_{j}. Because both *Ĥ* and $F^j$ are Hermitian as well as real operators, the correlation function [Eq. (2)] is a real and even function of *t*.

The main idea used in deriving the QI expression is to approximate the time integral in Eq. (1) using a steepest-descent (saddle-point) approximation taken around *t* = 0,

where the dot implies differentiation with respect to time. In order to use this approximation, it is necessary that $f\u0307$(0) = 0, which is guaranteed for even functions of *t*, and also that $f\u0308$(0) $<0$. This steepest-descent approximation makes it possible to express the full reaction rate in terms of expressions defined at *t* = 0, which can then be evaluated exactly with imaginary-time path-integral methods. The original QI expression [see Eq. (14) below] has been derived with one of two approaches: one which uses the energy-integral formulation of the reaction rate (with complicated transformations and two separate steepest-descent integrations)^{24} and another which involves multiplying and dividing by the delta-delta correlation function,

as the first step.^{29,40} Both derivations use the same prescription for choosing the dividing surfaces *x*_{1} and *x*_{2}.

Evaluating the trace in Eq. (2) in the position representation, we find the well-known expression^{1,40}

where “c.c.” denotes the complex conjugate, and

The delta-delta correlation function [Eq. (6)] is then simply

For brevity, we will drop the explicit dependence of the time correlation functions on the dividing surfaces in our notation. To apply the steepest-descent approximation [Eq. (4)] to Eq. (1), we set

which we note obeys the stationarity condition at time *t* = 0 because *Ċ*_{ff}(0) = 0. We are still free to choose the locations of the dividing surfaces, (*x*_{1}, *x*_{2}), and we should do this in such a way to ensure that *C*_{ff}(*t*) has a maximum at *t* = 0. Many choices for the dividing surfaces have been suggested, but some choices lead to a local minimum at zero time for *C*_{ff}(*t*), rather than a maximum.^{40} The original and most common choice is that the dividing surfaces obey

which typically leads to *C*_{ff}(*t*) being a maximum at *t* = 0. In our experience, there is always a solution to these equations with *x*_{1} = *x*_{2}, known as merged dividing surfaces, which describe either a minimum or first-order saddle point of *ρ*(*x*_{1}, *x*_{2}, 0). In the former case, another set of solutions exists as saddle points with *x*_{1} ≠ *x*_{2}, leading to what is known as split dividing surfaces. Adopting one of these choices of dividing surfaces such that *x*_{1} and *x*_{2} satisfy Eq. (10), we note that $Cff(0)=Cffsp(0)$, where

and thus suggest the approximation $Cff(t)\u2248Cffsp(t)$ for *t* ≠ 0.

We are now faced with a choice for performing the steepest-descent integration in time: we could either perform this approximate time integral on *C*_{dd}(*t*) and then take the spatial derivatives, or we could perform the integral directly on *C*_{ff}(*t*). Choosing first to perform the integral over time before evaluating the derivatives with respect to the positions of surfaces, the new function for the steepest-descent approximation and its second time derivative (in terms of the energy variance $\Delta Hdd2$) are

where *Ċ*_{dd}(0) = 0 and

The QI rate is then defined as

where we have assumed that the spatial derivatives of Δ*H*_{dd} can be neglected as they are much smaller than the derivatives of *C*_{dd}(0). Provided that the ratio *C*_{ff}(*t*)/*C*_{dd}(*t*) is slowly varying in time, Eq. (14) can be generalized for dividing surfaces that do not obey the condition in Eq. (10) such that^{40}

where we have slightly abused the notation and kept the QI subscript in this final expression.

An alternative approach is derived by employing a saddle-point approximation directly to the time integral of *C*_{ff}(*t*). In order to use this approximation, there is no particular requirement on the choice of dividing surfaces, provided *C*_{ff}(*t*) is a maximum at *t* = 0 for that choice. The idea of this direct integration has previously been introduced in the context of creating a “cumulant expansion,” where the expansion is truncated at the second order term^{41} (although it could be extended to higher order terms and, in addition, the exact high-temperature behavior could be incorporated analytically).^{41–43} The steepest-descent function is then

where

which defines the energy variance $\Delta Hff2$. The second-order cumulant expansion (2OCE) expression that results from this procedure is very similar to Eq. (15), namely,

We will show that both of these quantum instanton approximations have similar properties and break down for asymmetric systems for similar reasons.

### B. Behavior of the quantum instanton method for asymmetric systems

To test the approximations, it is useful to apply the QI and 2OCE methods to a simple 1D Eckart barrier,

where, choosing the same parameters as in Ref. 24, *V*_{0} = 0.425 eV, $a=1.36\u2009a0\u22121$, and the mass is 1060 *m*_{e} (*a*_{0} is the Bohr radius and *m*_{e} denotes the electron mass). *α* is the dimensionless asymmetry parameter controlling the degree of exothermicity between the reactants and products such that the barrier is symmetric for *α* = 1. We investigated a variety of different temperatures and asymmetries, and found that the QI method shows pathological behavior for low temperatures and large asymmetries. We will demonstrate this behavior throughout this paper, unless otherwise stated, with a particular choice of *T* = 100 K and *α* = 1.425, for which the QI prediction is an order of magnitude too large compared to the exact rate.

It would be reasonable to think that the QI method could be improved if a more suitable choice of dividing surfaces were found. We will therefore consider a range of choices for the dividing surfaces, as listed in Table I, with the above choice of parameters, showing examples of both split saddle points and a merged minimum of *C*_{dd}(0) and a merged saddle point of *C*_{ff}(0).^{44} The positions of the dividing surfaces are also indicated in Fig. 1, showing the behavior of the propagator and the flux-flux correlation function as functions of the dividing surface locations. It is possible to choose dividing surfaces that are not at special points,^{40} but we will show that it is not generally possible to find dividing surfaces that fix all the problems of the QI and 2OCE methods without needing to evaluate quantities at *t* ≠ 0.

Choice of surfaces . | Shorthand . | x_{1}
. | x_{2}
. |
---|---|---|---|

Split saddle points of C_{dd} | A | −0.968 | 0.063 |

Merged minimum of C_{dd} | B | −0.379 | −0.379 |

Merged saddle point of C_{ff} | C | −0.623 | −0.623 |

Semiclassical turning points | D | −1.847 | 0.828 |

Choice of surfaces . | Shorthand . | x_{1}
. | x_{2}
. |
---|---|---|---|

Split saddle points of C_{dd} | A | −0.968 | 0.063 |

Merged minimum of C_{dd} | B | −0.379 | −0.379 |

Merged saddle point of C_{ff} | C | −0.623 | −0.623 |

Semiclassical turning points | D | −1.847 | 0.828 |

To summarize Sec. II A, the QI and 2OCE methods essentially approximate the time dependence of the relevant time correlation functions as Gaussians. The functions *C*_{dd}(*t*) and *C*_{ff}(*t*) are shown in Fig. 2, where the surfaces are chosen to be the split surfaces of *ρ* if they exist or the merged surfaces of *ρ* otherwise. For increasing asymmetry *α*, it is clear that the functions become less and less Gaussian (a similar situation to that discussed in Refs. 45 and 46 for electron-transfer rates). This will lead to a significant error in the prediction of the rate constant. In fact, we find that none of the choices considered for the dividing surfaces make much difference in the QI prediction for the rate, as the QI prediction still deviates significantly from the exact value. In order to gain insight into this behavior, in Sec. III, we analyze the correlation functions in terms of semiclassical paths and discover the cause of the breakdown of the approximation.

## III. SEMICLASSICAL ANALYSIS

In this section, we will examine the various terms in the QI and 2OCE approximations using a semiclassical analysis, which gives the asymptotic behavior in the *ℏ* → 0 limit of a quantum-mechanical expression in terms of minimum-action paths. We will identify the dominant contributions that cause the observed deviation from the semiclassical instanton expression and determine what improvements would be necessary to fix these problems.

### A. Semiclassical contributions to quantum instanton quantities

All terms in the QI and 2OCE rate expressions are defined at *t* = 0, and thus, all of our analysis can be performed using the imaginary-time quantum propagator, *ρ*(*x*_{1}, *x*_{2}, 0), and its derivatives. Its semiclassical limit is the imaginary-time van Vleck propagator, which involves a sum over all minimum-action paths which travel from *x*_{1} to *x*_{2} in imaginary time *βℏ*/2. These are clearly the paths which dominate in a path-integral Monte Carlo evaluation of the quantities, in which paths are weighted by the exponential of the negative action. The minimum-action paths follow imaginary-time classical trajectories, which are equivalent to Newtonian trajectories on the upside-down potential-energy surface.^{47} In fact, there are two such minimum-action paths, one which bounces to the left and one to the right. Note that the imaginary-time classical trajectories which travel directly from one dividing surface to the other or which bounce more than once are saddle points of the action^{48} and do not contribute to the semiclassical limit. The semiclassical limit of the quantum propagator is thus given by

where $\u223c\u210f\u21920$ denotes the asymptotic behavior in the limit *ℏ* → 0,^{49} while *K*_{ℓ} and *K*_{r} are the semiclassical contributions from the left and right paths, given by

where *γ* ∈ {*ℓ*, *r*}. These propagators are functions of *x*_{1}, *x*_{2} and *τ*_{γ}, and in this case, both have the same imaginary-time length *τ*_{γ} = *βℏ*/2. The paths are written in terms of the classical Euclidean action (the action in imaginary time),

where *X*_{γ}(*τ*′) are the imaginary-time-dependent left and right paths from *x*_{1} to *x*_{2} via a “bounce.”^{50} Each action is related to an eikonal *W*_{γ} (also known as the reduced action) by the Legendre transform,^{17,48}

and

where *x* is integrated along the relevant path *X*_{γ}. The energy, *E*_{γ}, is chosen to solve $\u2202W\gamma \u2202E\gamma =\u2212\tau \gamma $ and the path bounces at the point where *E*_{γ} = *V* (*x*).

The two minimum-action paths are shown for the A, B, and D dividing surfaces in Fig. 3. It is important to notice that the two paths do not join together to form the semiclassical instanton solution, which is an imaginary-time periodic orbit with constant energy. Because the reaction is exothermic, for all choices of dividing surfaces, the left-bouncing path (green) follows half of the instanton periodic orbit, but the right-bouncing path (red) is spurious and has a much lower energy, *E*_{r}. This remains true even when using dividing surfaces placed at the turning points of the instanton trajectory (D), as advocated in the appendix of Ref. 24. In all these cases, the second half of the instanton periodic orbit cut at the location of the dividing surfaces is a first-order saddle point of the action, as it passes through a conjugate point,^{48} which is why the minimum-action path must follow a different classical trajectory. In fact, we found that for this system below about 142 K there is no way to split the instanton into two trajectories of equal imaginary time without encountering a conjugate point. The picture which we show is inconsistent with the qualitative picture suggested in the appendix of the original QI paper^{24} and, as we will show, the spurious path is the crux of the observed error of the QI method in asymmetric systems.

We use this semiclassical analysis in terms of minimum-action paths to find the dominant contributions to the various quantities used in the QI and 2OCE approximations. For instance, the semiclassical limit of *C*_{dd}(0) is simply

The semiclassical limits of the spatial first derivatives and mixed second derivative of the imaginary-time propagator are

where $pj\gamma =+2m[V(xj)\u2212E\gamma ]$ for *j* ∈ {1, 2}. Note that we take the positive root such that $pj\gamma $ is always a positive scalar and is therefore the magnitude of the momentum (without the direction). It is also clear from Fig. 3 that, contrary to what was previously thought,^{24,29} the dividing surfaces according to choice A in Table I are not close to the turning points of the semiclassical instanton, and under a semiclassical analysis, there is no relation between the two choices.

It is important to note that the terms proportional to $K\u21132$ and $Kr2$, which appear in *C*_{dd}(0), completely cancel out in the semiclassical limit of *C*_{ff}(0) as a consequence of the flux operators.

The real time derivatives (denoted by dots above quantities, as above) of the semiclassical propagator can be written using the Cauchy-Riemann equations, $K\u0307\gamma =i\u2202K\gamma \u2202\tau \gamma $ and $\u0116\gamma =i\u2202E\gamma \u2202\tau \gamma $, as

and

The combination of these two equations and Eq. (13) gives the energy variance of *C*_{dd} as

A similar analysis gives Δ*H*_{ff} in terms of the semiclassical quantities

The final expressions for the semiclassical limits of the QI and 2OCE rates are then simply obtained by inserting either Eqs. (27) and (30) into Eq. (15) for the QI result or Eqs. (27) and (31) into Eq. (18) for the 2OCE result.

Note that for less asymmetric systems or at higher temperatures (above 142 K in this case), it is possible to find a new definition for the dividing surface for which the two minimum-action paths have the same energy and thus describe the instanton orbit, as shown in Fig. 4. However, despite the fact that these two paths together describe the correct periodic orbit, we find that the QI approximation nonetheless overpredicts the rate by many orders of magnitude. This is because the left-bouncing path has a much smaller action than the right-bouncing path. As a result, the semiclassical analysis of *C*_{ff}(0) given in Eq. (27) is no longer a good estimate of the quantum-mechanical value. The breakdown of the semiclassical analysis is due to the fact that the term proportional to $K\u21132$ only cancels to first order in *ℏ* and cannot be neglected when compared with the much smaller *K*_{ℓ}*K*_{r} term. This problem is avoided by choosing dividing surfaces according to the standard prescriptions in Table I, for which *K*_{ℓ} and *K*_{r} are of a similar order of magnitude, for which the semiclassical analysis is valid.

For a completely symmetric system at any temperature, it is always possible to split the instanton in the middle, thereby generating two minimum-action paths of equal energy with *K*_{ℓ} = *K*_{r}. In the following, we shall use this fact to explain the success of the QI method for symmetric systems.

### B. Connection to the semiclassical instanton theory

The semiclassical limits to the QI and 2OCE approximations derived above show that the predicted rates do not just depend on the instanton path, but instead have contributions from the spurious path, *X*_{r}. However, the SCI result for the rate depends only on the instanton periodic orbit.^{17} The implication is that, contrary to the original conjecture, the SCI method is not a simple semiclassical approximation to the QI or 2OCE methods. In fact, it can be seen from both the above expressions and Fig. 3 that the contributions from the unphysical right-hand path will dominate in many situations and lead to erroneous rates. The only exception, when using the standard dividing surfaces satisfying Eq. (10), is for perfectly symmetric barriers where the two paths combine into the instanton periodic orbit.

These problems would go away if *ρ* were dominated by the semiclassical paths which together make up the instanton periodic orbit. This is of course already the case for a symmetric system. Here, we will show that if this were generally the case, then the SCI rate formula would be recovered as the asymptotic limit of the QI expression.

If the two paths together describe a periodic orbit, then we can use *E*_{ℓ} = *E*_{r} and hence $p1\u2113=p1r=p1$ and $p2\u2113=p2r=p2$. The flux-flux correlation function [Eq. (27)] thus simplifies to

and

In order to study the QI approximation based on Δ*H*_{dd}, we define dividing surfaces which obey Eq. (10) in the semiclassical limit [see also Eq. (26)], which implies that *K*_{ℓ} = *K*_{r}. The semiclassical limit of the energy variance Δ*H*_{dd} [Eq. (30)] becomes identical to that of Δ*H*_{ff} given in Eq. (33).

Equations (32) and (33) can thus be used to obtain the semiclassical limit of the QI and 2OCE approximations. Using the known expressions (for the 1D case)^{48}

and

Eq. (32) becomes

Using^{16}

where *E* is the energy of the instanton orbit, we obtain

which matches Miller’s original SCI expression^{14} in the one-dimensional case, bearing in mind that the total action of the instanton periodic orbit is given by *S*_{ℓ} + *S*_{r}. This equation is equal to the ImF instanton expression,^{18} which can be written in a number of equivalent ways^{17,19} including the ring-polymer formulation.^{21} Thus, if the right-hand and left-hand paths combine to form the instanton and the dividing surfaces are chosen to satisfy Eq. (10), the QI and 2OCE expressions reduce to the SCI rate in the semiclassical limit. In general, however, the two paths only join exactly into the instanton periodic orbit for a symmetric system; hence, the QI method is a particularly accurate method for symmetric barriers. For asymmetric barriers, the QI rate cannot be directly connected to the SCI expression.

## IV. PROJECTED QI

In this section, we suggest a modification of the quantum instanton approach to avoid sampling the spurious paths, which we will call the Projected Quantum Instanton (PQI) method. This defines a novel method in the spirit of the quantum instanton but which should give accurate rate predictions even for strongly asymmetric systems in the low-temperature limit.

In order to fix the QI method for asymmetric systems, it will be necessary to require that the two minimum-action paths have matching energies and thus combine into the instanton periodic orbit. Yet, for very asymmetric systems, it is impossible to choose dividing surfaces such that both minimum-action paths of imaginary-time *βℏ*/2 have the same energy. Therefore, we will need to relax the requirement that the left- and right-bouncing paths have equal imaginary-time lengths^{17} and will therefore also need a projection scheme which can categorize any general path into left and right sets. This can be achieved by defining the projected propagator

where the projection operators are written in terms of the Heaviside step function, *θ*, as

such that

Paths projected in this way can thus be categorized depending on whether their central point is to the right or left of the point *x*_{0}. The location of *x*_{0} is somewhat arbitrary for the following arguments as long as it appears between the instanton turning points. This way, ⟨*x*_{2}|*Û*_{ℓ}(−*iτ*_{ℓ})|*x*_{1}⟩ will be dominated by the left-bouncing minimum-action path and ⟨*x*_{2}|*Û*_{r}(−*iτ*_{r})|*x*_{1}⟩ by the right-bouncing minimum-action path. The projection operator we have used is the simplest choice to implement in our simulations. However, there are many other possible definitions which also pick out the correct semiclassical pathways. Note that the multidimensional extension of the approach follows directly by projecting along a reaction coordinate.

By inserting projection operators into the Miller–Schwartz–Tromp formula for the rate [Eq. (1)], we will effectively neglect contributions from pairs of paths which both bounce either to the left or to the right. This is expected to be a good approximation because within a semiclassical analysis these paths would give a zero contribution to the rate.^{16,17} We therefore propose to approximate the rate as

where

which is a general (nonsymmetric and complex) correlation function whose form permits the imaginary time to be different for the two propagators, although their sum, *τ*_{ℓ} + *τ*_{r} = *βℏ*, is fixed. This extra flexibility allows us to choose an appropriate dividing surface for which both the left- and right-bouncing paths can be minimum-action paths and also join together to describe the instanton. This idea follows a similar approach to that used in a first-principles derivation of the SCI method.^{16,17} It is absolutely necessary to project onto left and right paths when using a nonsymmetric split of the imaginary time to avoid finding the spurious minimum-action paths of a right-bouncing path in time *τ*_{ℓ} and vice versa. Note that the factor of 2 in Eq. (43) accounts for the alternative ordering of the projection operators, which integrates to the same result.^{17}

Equation (42) is an approximation in general, but remarkably it reproduces the exact rate of a free particle as shown in Appendix A. In this important limit, despite the fact that the projected and unprojected correlation functions are not equivalent, their integrals over time are identical.

The final expression for the PQI rate is the second-order cumulant expansion of

where the projected matrix elements are given by

and their time derivatives calculated using $\rho \u0307\gamma =i\u2202\rho \gamma \u2202\tau \gamma $ and $\rho \u0308\gamma =\u2212\u22022\rho \gamma \u2202\tau \gamma 2$. The PQI rate expression is then defined identically to the 2OCE approximation [Eq. (18)], except that *C*_{ff} is replaced by $CffP$ throughout.

By following a similar approach as in Sec. III B we obtain the semiclassical limit of $CffP(0)$ equal to that of Eq. (32), and of $\Delta Hff2$ given by Eq. (33). The only difference with the previous analysis is that in PQI the minimum-action paths really do describe the instanton periodic orbit and no extra assumptions need be made. One can therefore show that the semiclassical limit of the rate expression reduces to SCI using the results of Sec. III B. This result is actually unsurprising, as the derivation of SCI follows similar lines of reasoning in Ref. 17.

In principle, one can use any choice of dividing surfaces for which the two halves of the instanton periodic orbit are both minimum-action paths. We chose to place both the dividing surfaces as well as *x*_{0} at the location of the barrier maximum and found that this choice obeyed the rule in each case tested. One should then optimize *τ*_{ℓ} (keeping *τ*_{ℓ} + *τ*_{r} = *βℏ* fixed) until $\u010affP(0)=0$. However, the value of *τ*_{ℓ} obtained directly from SCI was found to be an excellent approximation to this optimal value.

We note that the steps used to derive the PQI method share some similarity to Wolynes’ nonadiabatic quantum-transition theory.^{51} In this approach, the two paths are forced to bounce left or right depending on a projection onto the two electronic states and for asymmetric systems may have different imaginary times. The idea for deriving new rate theories by ensuring that the instanton is the dominant path contributing to the rate has also been used in previous work on QTST^{21,35,36} and most recently in Ref. 52 for nonadiabatic rates.

The PQI method is also applicable to multidimensional problems and can be efficiently computed using the path-integral Monte Carlo approach with a simple extension to the standard methodology. The central bead of each path plays the role of the *x*′ variable in Eq. (40) and, for example, only paths for which *x*′ < *x*_{0} should contribute to *ρ*_{ℓ}. In this work, however, we implement the projection by integrating numerically over the allowed range of *x*′ and evaluate the imaginary-time propagator, as described in Appendix B.

## V. RESULTS AND DISCUSSION

In order to validate our semiclassical analysis of the quantum instanton and cumulant expansion, it is useful to compute the numerical values of the various terms making up these approximations and to compare the values obtained from quantum mechanics (using the eigenfunctions for the Eckart barrier discussed in Appendix B) and from the semiclassical approximation. The results of the quantum and semiclassical calculations for the same system and temperature as those used in Figs. 1 and 3 are compared in Table II. The semiclassical limit is seen to be very close to the quantum values for each of these choices of dividing surfaces, being at most a factor of two different. The discrepancy is of course larger for the dividing surface choice made in Fig. 4, which leads to *K*_{ℓ} ≫ *K*_{r}. The good agreement is also seen in Table III, where the QI and 2OCE rate predictions are compared using both quantum and semiclassical calculations. However, all these approximations are an order of magnitude larger than the exact rate constant, showing that the source of the error is in the QI or 2OCE rate formulas themselves. (The system and temperature were in fact specifically chosen to demonstrate this order of magnitude error.) This justifies our use of a semiclassical evaluation of the relevant quantities for analyzing the QI results. The PQI method predicts a rate of 4.34 × 10^{−12} a.u. without optimization of the imaginary-time split (i.e., using the semiclassical ratio of *τ*_{r}/*τ*_{ℓ} ≈ 0.258), but when it is optimized to *τ*_{r}/*τ*_{ℓ} ≈ 0.298, the rate prediction is 4.85 × 10^{−12} a.u., which is even closer to the exact result (cf. the caption of Table III).

. | C_{dd}(0)
. | C_{ff}(0)
. | ΔH_{dd}
. | ΔH_{ff}
. | ||||
---|---|---|---|---|---|---|---|---|

Div. surf. . | QM . | SC . | QM . | SC . | QM . | SC . | QM . | SC . |

A | 5.40 × 10^{−9} | 4.16 × 10^{−9} | 4.60 × 10^{−14} | 3.88 × 10^{−14} | 0.003 23 | 0.003 23 | 0.004 37 | 0.004 48 |

B | 5.12 × 10^{−9} | 4.03 × 10^{−9} | 6.60 × 10^{−14} | 5.59 × 10^{−14} | 0.003 37 | 0.003 34 | 0.004 54 | 0.004 60 |

C | 1.44 × 10^{−7} | 1.31 × 10^{−7} | 4.11 × 10^{−14} | 3.07 × 10^{−14} | 0.001 29 | 0.001 06 | 0.004 09 | 0.004 41 |

D | 1.58 × 10^{−8} | 2.27 × 10^{−8} | 1.61 × 10^{−14} | 3.00 × 10^{−14} | 0.002 31 | 0.002 95 | 0.003 51 | 0.004 03 |

. | C_{dd}(0)
. | C_{ff}(0)
. | ΔH_{dd}
. | ΔH_{ff}
. | ||||
---|---|---|---|---|---|---|---|---|

Div. surf. . | QM . | SC . | QM . | SC . | QM . | SC . | QM . | SC . |

A | 5.40 × 10^{−9} | 4.16 × 10^{−9} | 4.60 × 10^{−14} | 3.88 × 10^{−14} | 0.003 23 | 0.003 23 | 0.004 37 | 0.004 48 |

B | 5.12 × 10^{−9} | 4.03 × 10^{−9} | 6.60 × 10^{−14} | 5.59 × 10^{−14} | 0.003 37 | 0.003 34 | 0.004 54 | 0.004 60 |

C | 1.44 × 10^{−7} | 1.31 × 10^{−7} | 4.11 × 10^{−14} | 3.07 × 10^{−14} | 0.001 29 | 0.001 06 | 0.004 09 | 0.004 41 |

D | 1.58 × 10^{−8} | 2.27 × 10^{−8} | 1.61 × 10^{−14} | 3.00 × 10^{−14} | 0.002 31 | 0.002 95 | 0.003 51 | 0.004 03 |

. | k_{QI} × 10^{12}
. | k_{2OCE} × 10^{12}
. | ||
---|---|---|---|---|

Div. surf. . | QM . | SC . | QM . | SC . |

A | 54.6 | 46.0 | 40.3 | 33.2 |

B | 75.1 | 64.1 | 55.7 | 46.6 |

C | 123 | 111 | 38.6 | 26.7 |

D | 26.7 | 39.0 | 17.5 | 28.5 |

. | k_{QI} × 10^{12}
. | k_{2OCE} × 10^{12}
. | ||
---|---|---|---|---|

Div. surf. . | QM . | SC . | QM . | SC . |

A | 54.6 | 46.0 | 40.3 | 33.2 |

B | 75.1 | 64.1 | 55.7 | 46.6 |

C | 123 | 111 | 38.6 | 26.7 |

D | 26.7 | 39.0 | 17.5 | 28.5 |

The example system used in Fig. 3 shows a significant error only at low enough temperatures. However, we stress that the error seen in the QI method depends on both the asymmetry and the temperature such that for more asymmetric systems, even high temperatures can result in very large errors. For a more general outlook, Fig. 5 shows the relative error between the exact rate and the rates calculated using the four methods (QI, 2OCE, SCI, and PQI) as a function of asymmetry, *α*, and for higher temperatures than that considered in Tables II and III.^{53} To generate these results, we have generally favored A-type dividing surfaces, although it would have been possible to choose surfaces for the 2OCE calculation that minimize the contributions of higher order terms.^{41} Nonetheless, our semiclassical analysis implies that no dividing surface choice will significantly improve these results.

The temperature dependence of the exact and approximate rate constants is shown in Fig. 6(a), along with the corresponding relative errors in Fig. 6(b). While the SCI rates tend to a small constant relative error in the deep tunneling regime,^{17,54} the error in the QI rates increases exponentially with asymmetry. We attribute this growing error in the QI rate to the increasing contribution of the right-hand path in Fig. 3. In principle, 2OCE could be improved by going to a higher order expansion;^{41} however, in practice, very high orders will typically be required to cancel the spurious dominant contribution at zero time. In contrast, the PQI results in Figs. 5 and 6 show a substantial improvement over both the QI and 2OCE results and are at least as accurate as SCI. For these PQI results, the value of *τ*_{ℓ} is not optimized but simply chosen to match the imaginary time for the left bounce in the instanton trajectory.

Close to the crossover temperature, *T*_{c} = *ℏω*_{b}/2*πk*_{B} (where *ω*_{b} is the absolute value of the imaginary barrier frequency),^{17} the SCI approximation breaks down. The accuracy of the PQI rate is also decreased at higher temperatures, although it is not much worse than the original QI and 2OCE methods, a statement which we can quantify by studying the rate of a free particle (and hence the high-temperature limit for arbitrary barriers). In Appendix A, we calculate the PQI rate for this case and our analytical results explain why the PQI prediction slightly underestimates the exact rate in the high-temperature regime. However, we also show that direct integration of $CffP(t)$ yields the exact rate for the free particle, suggesting that the only error lies in the steepest-descent approximation.

Overall, these promising results show that the new PQI method is a substantial improvement over the traditional QI and 2OCE methods. Future extensions of PQI for use with path-integral methods are therefore likely to fulfill the initial promise of QI as an extension to SCI, suitable for accurately predicting deep-tunneling reaction rates in highly anharmonic and asymmetric systems.

## VI. CONCLUSIONS

In conclusion, we have presented a semiclassical analysis of the quantum instanton method. We have shown that the dominant contribution to the QI expression can arise from spurious paths, which leads to very large errors for very asymmetric barriers, especially at low temperatures. Consequently, despite conventional wisdom suggesting otherwise, no choice of dividing surface can remove this problem. We justify our analysis by showing that semiclassical evaluation of the relevant quantities yields very similar numerical values to those from exact quantum mechanics for the asymmetric Eckart barrier. The major discrepancy from the exact rate therefore lies with the QI and 2OCE approximations themselves.

From our analysis, we conclude that the SCI rate is the semiclassical limit of the QI rate *only* when the barrier is perfectly symmetric. Therefore, the QI method provides an accurate prediction of the quantum rate for symmetric barriers. However, we find that SCI is much more accurate than QI at describing tunneling in asymmetric systems, although the SCI clearly cannot describe anharmonicity as well as the QI does.

Our proposed new method, PQI, fixes the inherent path-sampling problems with the QI and 2OCE methods mentioned above, and, at the same time, takes into account anharmonicities that SCI neglects. There exists some similarity between the PQI and theories employed to simulate electron-transfer reactions,^{51,52} and we hope that this study will inspire future development of nonadiabatic rate theories. It remains to be seen whether a path-integral application of PQI will become a competitive method for reaction rate calculations of complex systems. Although the theory can be easily generalized to multidimensional systems, it may not be so easy to apply to atomistic simulations of reactions in solution where SCI theory is not well defined.^{15} RPMD rate theory^{7} also dominantly samples instanton paths even for asymmetric systems^{21} and, in addition, includes classical recrossing effects neglected by SCI, QI, and PQI. It thus remains the method of choice for atomistic simulations of reactions in solution.^{55}

## ACKNOWLEDGMENTS

The authors acknowledge support from the Swiss National Science Foundation through the NCCR MUST (Molecular Ultrafast Science and Technology) Network and Dr. Konstantin Karandashev and Joseph Lawrence for useful discussions. M.J.T. is supported by an ETH Zurich Research Grant.

### APPENDIX A: PQI, 2OCE, AND QI RATES FOR THE FREE PARTICLE

It is informative to examine the behavior of PQI for the case of the free particle, which can be treated analytically. The exact propagator for the free particle is well known,^{56}

The integrals over *x*′ can be done analytically with standard results^{57} to obtain

where Φ(*z*) is the Gauss error function and

Note that, as required from Eq. (41), *ρ* = *ρ*_{ℓ} + *ρ*_{r}.

The free-particle delta-delta correlation function is

and the stationarity condition [Eq. (10)] for the dividing surfaces, equivalent to setting *∂C*_{dd}(0)/*∂x*_{γ} = 0, implies *x*_{1} = *x*_{2}. Taking *x*_{0} = *x*_{1} = *x*_{2}, we find *ρ*_{ℓ} = *ρ*_{r} = *ρ*/2.

For *x*_{1} = *x*_{2}, the exact quantum flux-flux correlation function for the free particle is^{1}

Quantities required for the QI rate [Eq. (14)] are thus $\Delta Hdd=2/\beta $ and $Cff(0)=(\pi \u210f2\beta 2)\u22121$. The resulting QI rate,

is a factor of $\pi /2$ from the correct free-particle result, *kQ*_{r} = (2*πℏβ*)^{−1}, leading to an overestimation of the rate by 25%.^{24}

For the 2OCE rate (18), the relevant quantity is instead $\Delta Hff=\u20096/\beta $, giving the rate

which differs from the exact free particle rate by a factor $\pi /6$ (an underestimation by about 28%).

Likewise, we can calculate the projected correlation function for the free particle using Eqs. (44) and (A3) with the substitution *τ*_{γ} → *τ*_{γ} + *it* to obtain

The relevant quantities needed to calculate the PQI rate for the free particle are therefore

and

Therefore, the PQI rate for the free particle is given as

which underestimates the exact rate by about 37%.

It is possible to directly integrate $CffP(t)$ over *t* to yield the exact rate for this system. Then, using Eq. (42), we get *kQ*_{r} = (2*πℏβ*)^{−1}, which is the exact rate for this system.^{24} Therefore, for the special case of the free particle, no approximation is made by adding projections in Eq. (42). However, in their simplest forms, QI, 2OCE, and PQI each give a different result for the rate of a free particle, none of which is correct. The errors in these approximate methods stem from the fact that neither the standard nor the projected flux-flux correlation function is a Gaussian function of *t*, and thus, the steepest-descent approximation introduces an error. This is well understood from the original QI work, and simple corrections have been suggested^{24,41,58} which could also be applied to improve the PQI rate formula.

### APPENDIX B: ECKART BARRIER EIGENFUNCTIONS

In this appendix, we derive closed-form expressions for the wavefunctions of the one-dimensional Eckart barrier. Eckart^{59} presented a derivation for only one wavefunction, which was enough to obtain the reflection and transmission coefficients. For our purpose, however, we will need two orthogonal wavefunctions to serve as a complete basis. In nearly all other respects, we follow his approach.

The Schrödinger equation for this problem can be written as

where comparison with Eq. (19) gives *ξ* = −*e*^{2πx/l}, *A* = *V*_{0}(1 − *α*), $B=V0(1+\alpha )2$, and *l* = *π*/*a*. Following Eckart, we define

and then rearrange into the form

Note that if *E* < 0, the root in Eq. (B2a) is taken such that $Im\u2009\alpha \xaf<0$. We assume that *A* < 0, i.e., the asymptote on the right-hand side is lower than that on the left.

It is known that the hypergeometric function solves Eq. (B3). Klein^{60} showed that solutions to this differential equation can be obtained using the ansatz,

with

However, *F*[*a*_{0}, *b*_{0}, *c*_{0}, *ξ*] is only one possible solution to the hypergeometric differential equation, and instead we choose to use Forsyth’s solutions^{61} called III and XIII to give

The former is only a solution if *E* > *A* and the second if *E* > 0. In order to find the *ξ* → −*∞* limit, we first apply Abramowitz and Stegun’s formula 15.3.7,^{62} and after premultiplying by a constant, we find the desired asymptotic limits for *u*_{1} to be

where

and $k\alpha \xaf=2\pi \alpha \xaf/l$ and $k\beta \xaf=2\pi \beta \xaf/l$. The asymptotic limits of the second wavefunction are (using Abramowitz and Stegun 15.3.9)^{62}

where

Note that Eckart’s results for reflection and transmission are given by |*B*_{2}/*A*_{2}|^{2} and $|\beta \xaf/\alpha \xaf|/|A2|2$ given a particle incident on the left.

Finally, we normalize the wavefunctions

such that

Using

we find the normalization constants should be chosen as

where $k\alpha \xaf\u2032=\u2202k\alpha \xaf\u2202E=k\alpha \xaf/(2E)$ and $k\beta \xaf\u2032=\u2202k\beta \xaf\u2202E=k\beta \xaf/(2(E\u2212A))$. The two wavefunctions must be orthogonal because one corresponds to an incoming particle from the left and the other from the right. Mathematically, it occurs because $B2/k\alpha \xaf\u2032+C1*/k\beta \xaf\u2032=0$.

The propagator can be evaluated numerically as an integral over energies using

where *t* can also be complex in order to obtain the density matrix. Derivatives of the propagator can be evaluated by explicit differentiation of the wavefunctions, although in practice we found that it was simpler and sufficiently accurate to evaluate them with finite differences.

## REFERENCES

We were unable to locate an asymmetric system with split saddle points of *C*_{ff}(0).

The error in the QI rates is slightly larger than those reported in Ref. 25 because they used an *ad hoc* correction to the energy variance for their calculations, but this does not affect our conclusions.