Semiclassical analysis of the quantum instanton approximation

TL;DR

This paper analyzes the relationship between quantum and semiclassical instanton approximations for reaction rates, revealing limitations in asymmetric systems and proposing a modified approach for improved accuracy.

Contribution

It demonstrates the connection between quantum and semiclassical instanton methods and introduces a modification to improve predictions for asymmetric barriers.

Findings

Two dominant minimum-action paths in quantum instanton expression

Semiclassical instanton theory is accurate mainly for symmetric systems

Proposed modification avoids unphysical path sampling, improving accuracy

Abstract

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 only be obtained directly from the semiclassical limit of the quantum instanton 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.

Tables

Table 1. Table 1: Definition and positions of dividing surfaces for the asymmetric Eckart barrier with α = 1.425 𝛼 1.425 \alpha=1.425 and T = 100 𝑇 100 T=100 K. The remaining parameters are defined in the main text.

Choice of surfaces	Shorthand	$x_1$	$x_2$
Split saddle points of $C_{\mathrm{dd}}$	A	$-0.968$	$0.063$
Merged minimum of $C_{\mathrm{dd}}$	B	$-0.379$	$-0.379$
Merged saddle point of $C_{\mathrm{ff}}$	C	$-0.623$	$-0.623$
Semiclassical turning points	D	$-1.847$	$0.828$

Table 2. Table 2: Numerical values (in a.u.) of various quantities used in the QI expressions are calculated using quantum-mechanical (QM) and semiclassical (SC) methods for the specified dividing surfaces. The parameters of the potential are the same as in Table 1 and the temperature is 100 K. In each case, the semiclassical version is within a factor of two of the quantum result, confirming that our semiclassical analysis of the quantum instanton is valid.

Div. Surf.	$C_{\mathrm{dd}}\left(0\right)$		$C_{\mathrm{ff}}\left(0\right)$		$\Delta H_{\mathrm{dd}}$		$\Delta H_{\mathrm{ff}}$
	QM	SC	QM	SC	QM	SC	QM	SC
A	$5.40\times {10}^{-9}$	$4.16\times {10}^{-9}$	$4.60\times {10}^{-14}$	$3.88\times {10}^{-14}$	0.00323	0.00323	0.00437	0.00448
B	$5.12\times {10}^{-9}$	$4.03\times {10}^{-9}$	$6.60\times {10}^{-14}$	$5.59\times {10}^{-14}$	0.00337	0.00334	0.00454	0.00460
C	$1.44\times {10}^{-7}$	$1.31\times {10}^{-7}$	$4.11\times {10}^{-14}$	$3.07\times {10}^{-14}$	0.00129	0.00106	0.00409	0.00441
D	$1.58\times {10}^{-8}$	$2.27\times {10}^{-8}$	$1.61\times {10}^{-14}$	$3.00\times {10}^{-14}$	0.00231	0.00295	0.00351	0.00403

Table 3. Table 3: Prediction of rates (in a.u.) using different methods, where k 2 ​ O ​ C ​ E subscript 𝑘 2 O C E k_{\mathrm{2OCE}} is the rate from the cumulant expansion, and k QI subscript 𝑘 QI k_{\mathrm{QI}} is the rate from the QI method. The parameters of the potential are the same as in Table 1 and the temperature is 100 K. The exact rate is 4.63 × 10 − 12 4.63 superscript 10 12 4.63\times 10^{-12} a.u. and the semiclassical instanton rate is 3.93 × 10 − 12 3.93 superscript 10 12 3.93\times 10^{-12} a.u. For comparison, the classical TST rate is 8.31 × 10 − 26 8.31 superscript 10 26 8.31\times 10^{-26} a.u., which confirms that quantum tunneling effects are of significant importance. Rate constants were evaluated using either the quantum-mechanical (QM) or semiclassical (SC) components from Table  2 .

Div. Surf.	$k_{\mathrm{QI}}\times {10}^{12}$		$k_{2OCE}\times {10}^{12}$
	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

Equations

$$k{Q}_{\mathrm{r}}=\frac{1}{2}{\int }_{-\infty }^{\infty }dt{C}_{\mathrm{ff}}(t).$$

$$C_{\mathrm{ff}}(t)=\mathrm{Tr}\left(e^{-\beta \hat{H}/2}{\hat{F}}_1{e}^{-\beta \hat{H}/2}{e}^{i\hat{H}t/\hslash }{\hat{F}}_2{e}^{-i\hat{H}t/\hslash }\right),$$

$${\hat{F}}_j=\frac{1}{2m}\left({\hat{\delta }}_j\hat{p}+\hat{p}{\hat{\delta }}_j\right)(j\in \{1,2\}),$$

$${\int }_{-\infty }^{\infty }dt{e}^{f(t)}\approx {\int }_{-\infty }^{\infty }dt{e}^{f(0)+\ddot{f}(0)t^2/2}=\sqrt{\frac{2\pi }{-\ddot{f}(0)}}{e}^{f(0)},$$

$$C_{\mathrm{dd}}(t)=\mathrm{Tr}\left(e^{-\beta \hat{H}/2}{\hat{\delta }}_1{e}^{-\beta \hat{H}/2}{e}^{i\hat{H}t/\hslash }{\hat{\delta }}_2{e}^{-i\hat{H}t/\hslash }\right),$$

$$C_{\mathrm{ff}}(t)={\left(\frac{\hslash }{2m}\right)}^2\left(\frac{{\partial }^2\rho }{\partial x_1\partial x_2}{\rho }^{\ast }-\frac{\partial \rho }{\partial x_1}\frac{\partial {\rho }^{\ast }}{\partial x_2}\right)+c.c.,$$

$$\rho \equiv \rho (x_1,x_2,t)=⟨x_2|e^{-\beta \hat{H}/2-i\hat{H}t/\hslash }|x_1⟩.$$

$$C_{\mathrm{dd}}(t)=|\rho (x_1,x_2,t){|}^2.$$

$$f_{\mathrm{ff}}(t)=\ln C_{\mathrm{ff}}(t),$$

$$\frac{\partial \rho (x_1,x_2,0)}{\partial x_1}=\frac{\partial \rho (x_1,x_2,0)}{\partial x_2}=0,$$

$$C_{\mathrm{ff}}^{\mathrm{sp}}(t)={\left(\frac{\hslash }{2m}\right)}^2\frac{{\partial }^2}{\partial x_1\partial x_2}{C}_{\mathrm{dd}}(t),$$

$$f_{\mathrm{dd}}(t)$$

$${\frac{d^2{f}_{\mathrm{dd}}}{d{t}^2}|}_{t=0}$$

$${\ddot{C}}_{\mathrm{dd}}(t)=2|\dot{\rho }{|}^2+\ddot{\rho }{\rho }^{\ast }+\rho {\ddot{\rho }}^{\ast }.$$

$$\begin{array}{rl}{k}_{\mathrm{QI}}{Q}_{\mathrm{r}}& =\frac{1}{2}{\left(\frac{\hslash }{2m}\right)}^2\frac{{\partial }^2{C}_{\mathrm{dd}}(0)}{\partial x_1\partial x_2}{\int }_{-\infty }^{\infty }dt{e}^{-\Delta H_{\mathrm{dd}}^2{t}^2/{\hslash }^2}\\ & =\frac{\hslash \sqrt{\pi }}{2\Delta H_{\mathrm{dd}}}{C}_{\mathrm{ff}}^{\mathrm{sp}}(0),\end{array}$$

$$k_{\mathrm{QI}}{Q}_{\mathrm{r}}=\frac{\hslash \sqrt{\pi }}{2\Delta H_{\mathrm{dd}}}{C}_{\mathrm{ff}}(0),$$

$$f_{\mathrm{ff}}(t)$$

$${\frac{d^2{f}_{\mathrm{ff}}}{d{t}^2}|}_{t=0}$$

$$\begin{array}{rl}{\ddot{C}}_{\mathrm{ff}}(t)=& {\left(\frac{\hslash }{2m}\right)}^2(\frac{{\partial }^2\rho }{\partial x_1\partial x_2}{\ddot{\rho }}^{\ast }+2\frac{{\partial }^2\dot{\rho }}{\partial x_1\partial x_2}{\dot{\rho }}^{\ast }+\frac{{\partial }^2\ddot{\rho }}{\partial x_1\partial x_2}{\rho }^{\ast }\\ & -\frac{\partial \ddot{\rho }}{\partial x_1}\frac{\partial {\rho }^{\ast }}{\partial x_2}-2\frac{\partial \dot{\rho }}{\partial x_1}\frac{\partial {\dot{\rho }}^{\ast }}{\partial x_2}-\frac{\partial \rho }{\partial x_1}\frac{\partial {\ddot{\rho }}^{\ast }}{\partial x_2})+c.c.,\end{array}$$

$$\begin{array}{rl}{k}_{2OCE}{Q}_{\mathrm{r}}& =\frac{1}{2}{C}_{\mathrm{ff}}(0){\int }_{-\infty }^{\infty }dt{e}^{-\Delta H_{\mathrm{ff}}^2{t}^2/{\hslash }^2}\\ & =\frac{\hslash \sqrt{\pi }}{2\Delta H_{\mathrm{ff}}}{C}_{\mathrm{ff}}(0).\end{array}$$

$$V(x)=\frac{(1-\alpha )V_0}{1+e^{-2ax}}+\frac{(1+\sqrt{\alpha }{)}^2{V}_0}{4{\cosh }^2(ax)},$$

$$\rho =⟨x_2⟩e^{-\beta \hat{H}/2}{x}_1\underset{\hslash \to 0}{\sim }{K}_{\ell }+K_r,$$

$$K_{\gamma }=(2\pi \hslash {)}^{-\frac{1}{2}}{\left(-\frac{{\partial }^2{S}_{\gamma }}{\partial x_1\partial x_2}\right)}^{\frac{1}{2}}{e}^{-S_{\gamma }/\hslash },$$

$$S_{\gamma }={\int }_0^{{\tau }_{\gamma }}d{\tau }^{\prime }\left\{\frac{1}{2}m{\left|\frac{d{X}_{\gamma }({\tau }^{\prime })}{d{\tau }^{\prime }}\right|}^2+V[X_{\gamma }({\tau }^{\prime })]\right\},$$

$$S_{\gamma }\equiv S_{\gamma }(x_1,x_2,{\tau }_{\gamma })=W_{\gamma }(x_1,x_2,E_{\gamma })+E_{\gamma }{\tau }_{\gamma },$$

$$W_{\gamma }\equiv W_{\gamma }(x_1,x_2,E_{\gamma })={\int }_{X_{\gamma }}dx\sqrt{2m[V(x)-E_{\gamma }]},$$

$$C_{\mathrm{dd}}(0)\underset{\hslash \to 0}{\sim }(K_{\ell }+K_r{)}^2=K_{\ell }^2+2K_{\ell }{K}_r+K_r^2.$$

$$\frac{\partial \rho }{\partial x_1}$$

$$\frac{\partial \rho }{\partial x_2}$$

$$\frac{{\partial }^2\rho }{\partial x_1\partial x_2}$$

$$C_{\mathrm{ff}}(0)\underset{\hslash \to 0}{\sim }\frac{1}{2m^2}(p_1^{\ell }+p_1^r)(p_2^{\ell }+p_2^r)K_{\ell }{K}_r.$$

$$\dot{\rho }\underset{\hslash \to 0}{\sim }-\frac{i}{\hslash }\left(E_{\ell }{K}_{\ell }+E_r{K}_r\right)$$

Full text

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 favoured.

A rigorous formulation of reaction rate theory is provided by the flux-side or flux-flux time correlation function formalism.Miller, Schwartz, and Tromp (1983) These correlation functions can be evaluated with approximate quantum-dynamics methods, including linearized semiclassical initial-value representation (LSC-IVR),Miller (2001a) centroid molecular dynamics,Cao and Voth (1994); Jang and Voth (1999) ring-polymer molecular dynamics (RPMD),Craig and Manolopoulos (2004, 2005a, 2005b) and the recently-developed Matsubara dynamicsHele et al. (2015a, b) and its approximations.Willatt, Ceriotti, and Althorpe (2018); Trenins and Althorpe (2018) 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.Miller (2001b) 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.Miller (1974); Pollak and Liao (1998)

Semiclassical instanton (SCI) theory provides one of the simplest formulations of a QTST for multidimensional tunneling.Miller (1975); Richardson (2018a) It can be derived rigorously by taking asymptotic ($\hbar\rightarrow 0$) approximations to the flux-flux correlation functionRichardson (2016a) and is known in a number of formulations,Richardson (2018b) including one obtained from the “ImF” approximation,Affleck (1981) which can all be shown to be equivalent.Althorpe (2011) 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 approachAndersson et al. (2009); Richardson and Althorpe (2009); Rommel, Goumans, and Kästner (2011) provides a computationally efficient algorithm for applying the theory in practice, and can be automated for application to complex systems.Vaillant, Althorpe, and Wales (2019)

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),Miller et al. (2003a) 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.Yamamoto and Miller (2004) The QI method has been used to calculate reaction rates in many multidimensional systems,Zhao, Yamamoto, and Miller (2004); Wang and Zhao (2009); Yamamoto and Miller (2005) and is a particularly efficient approach for the direct calculation of kinetic isotope effects.Vaníček et al. (2005); Karandashev et al. (2017); Karandashev and Vaníček (2015, 2017) 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.Karandashev (2017); Richardson (2018b)

There are also a number of other path-integral-based QTSTs developed to go beyond the SCI approximation. Unlike in classical mechanics,Chandler (1978) 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 theoryMills et al. (1997) and ring-polymer transition-state theory, which is itself related to the RPMD rate theory.Richardson and Althorpe (2009) 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.Hele and Althorpe (2013); Althorpe and Hele (2013) This provides an extension of the centroid-based QTST of Voth, Chandler and Miller,Voth, Chandler, and Miller (1989) a method which does not dominantly sample the instanton and thus fails for asymmetric barriers.Makarov and Topaler (1995)

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, that 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.

The 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 the following section 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.Miller et al. (2003a); Yamamoto and Miller (2004); Zhao, Yamamoto, and Miller (2004); Vaníček et al. (2005) Our conclusions remain valid for these multidimensional cases.

II.1 Derivation

We begin by rederiving the QI approximation from the Miller–Schwartz–Tromp expression for the reaction rate constantMiller, Schwartz, and Tromp (1983)

[TABLE]

Here, $Q_{\mathrm{r}}$ is the reactant partition function and $C_{\mathrm{ff}}$ is the symmetrized flux-flux time correlation function,

[TABLE]

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

[TABLE]

where $m$ is the mass and $\hat{\delta}_{j}=\delta(\hat{x}-x_{j})$ indicates a Dirac delta function centred at the $j$-th dividing surface located at $x=x_{j}$. Because both $\hat{H}$ and $\hat{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$,

[TABLE]

where the dot implies differentiation with respect to time. In order to use this approximation, it is necessary that $\dot{f}(0)=0$, which is guaranteed for even functions of $t$, and also that $\ddot{f}(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)Miller et al. (2003a) and another which involves multiplying and dividing by the delta-delta correlation function,

[TABLE]

as the first step.Vaníček et al. (2005); Aieta and Ceotto (2017) 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 expressionMiller, Schwartz, and Tromp (1983); Aieta and Ceotto (2017)

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

which we note obeys the stationarity condition at time $t=0$ because $\dot{C}_{\mathrm{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_{\mathrm{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_{\mathrm{ff}}(t)$, rather than a maximum.Aieta and Ceotto (2017) The original and most common choice is that the dividing surfaces obey

[TABLE]

which typically leads to $C_{\mathrm{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 $\rho(x_{1},x_{2},0)$. In the former case, another set of solutions exists as saddle points with $x_{1}\neq 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 $C_{\mathrm{ff}}(0)=C^{\mathrm{sp}}_{\mathrm{ff}}(0)$, where

[TABLE]

and thus suggest the approximation $C_{\mathrm{ff}}(t)\approx C^{\mathrm{sp}}_{\mathrm{ff}}(t)$ for $t\neq 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_{\mathrm{dd}}(t)$ and then take the spatial derivatives, or we could perform the integral directly on $C_{\mathrm{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 H_{\mathrm{dd}}^{2}$) are

[TABLE]

where $\dot{C}_{\mathrm{dd}}(0)=0$ and

[TABLE]

The QI rate is then defined as

[TABLE]

where we have assumed that the spatial derivatives of $\Delta H_{\textrm{dd}}$ can be neglected as they are much smaller than the derivatives of $C_{\mathrm{dd}}(0)$. Provided that the ratio $C_{\mathrm{ff}}(t)/C_{\mathrm{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 thatAieta and Ceotto (2017)

[TABLE]

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_{\mathrm{ff}}(t)$. In order to use this approximation, there is no particular requirement on the choice of dividing surfaces, provided $C_{\mathrm{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 Ceotto, Yang, and Miller (2005) (although it could be extended to higher order terms and, in addition, the exact high-temperature behavior could be incorporated analytically).Ceotto, Yang, and Miller (2005); Yang, Yamamoto, and Miller (2006); Predescu (2004) The steepest-descent function is then

[TABLE]

where

[TABLE]

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

[TABLE]

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

II.2 Behaviour 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,

[TABLE]

where, choosing the same parameters as in Ref. Miller et al., 2003a, $V_{0}=0.425$ eV, $a=1.36\;a_{0}^{-1}$, the mass is 1060 $m_{\mathrm{e}}$ ($a_{0}$ is the Bohr radius and $m_{\mathrm{e}}$ denotes the electron mass). $\alpha$ is the dimensionless asymmetry parameter controlling the degree of exothermicity between the reactants and products such that the barrier is symmetric for $\alpha=1$. We investigated a variety of different temperatures and asymmetries, and found that the QI method shows pathological behaviour for low temperatures and large asymmetries. We will demonstrate this behaviour throughout the paper, unless otherwise stated, with a particular choice of $T=100$ K and $\alpha=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, listed in Table 1, with the above choice of parameters, showing examples of both split saddle points and a merged minimum of $C_{\mathrm{dd}}(0)$ and a merged saddle point of $C_{\mathrm{ff}}(0)$. 111We were unable to locate an asymmetric system with split saddle points of $C_{\mathrm{ff}}(0)$. 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,Aieta and Ceotto (2017) 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\neq 0$.

To summarize the previous section, the QI and 2OCE methods essentially approximate the time dependence of the relevant time correlation functions as Gaussians. The functions $C_{\mathrm{dd}}(t)$ and $C_{\mathrm{ff}}(t)$ are shown in Fig. 2, where the surfaces are chosen to be the split surfaces of $\rho$ if they exist, or the merged surfaces of $\rho$ otherwise. For increasing asymmetry $\alpha$ it is clear that the functions become less and less Gaussian (a similar situation to that discussed in Refs. Huo, Miller III, and Coker, 2013; Richardson and Thoss, 2014 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 behaviour, in the next section 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 behaviour in the $\hbar\rightarrow 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.

III.1 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, $\rho(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 $\beta\hbar/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 minus the action. The minimum-action paths follow imaginary-time classical trajectories, which are equivalent to Newtonian trajectories on the upside-down potential-energy surface.Miller (1971) 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 Gutzwiller (1990) and do not contribute to the semiclassical limit. The semiclassical limit of the quantum propagator is thus given by

[TABLE]

where $\underset{\hbar\rightarrow 0}{\thicksim}$ denotes the asymptotic behaviour in the limit $\hbar\rightarrow 0$,Bender and Orszag (1978) while $K_{\ell}$ and $K_{r}$ are the semiclassical contributions from the left and right paths, given by

[TABLE]

where $\gamma\in\{\ell,r\}$. These propagators are functions of $x_{1},x_{2}$ and $\tau_{\gamma}$, and in this case both have the same imaginary-time length $\tau_{\gamma}=\beta\hbar/2$. The paths are written in terms of the classical Euclidean action (the action in imaginary time)

[TABLE]

where $X_{\gamma}(\tau^{\prime})$ are the imaginary-time-dependent left and right paths from $x_{1}$ to $x_{2}$ via a “bounce”.Richardson, Bauer, and Thoss (2015) Each action is related to an eikonal $W_{\gamma}$ (also known as the reduced action) by the Legendre transform Gutzwiller (1990); Richardson (2018b)

[TABLE]

and

[TABLE]

where $x$ is integrated along the relevant path $X_{\gamma}$. The energy, $E_{\gamma}$, is chosen to solve $\frac{\partial{W_{\gamma}}}{\partial{E_{\gamma}}}=-\tau_{\gamma}$ and the path bounces at the point where $E_{\gamma}=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. Miller et al., 2003a. 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,Gutzwiller (1990) 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 paperMiller et al. (2003a) 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_{\mathrm{dd}}(0)$ is simply

[TABLE]

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

[TABLE]

where $p_{j}^{\gamma}=+\sqrt{2m[V(x_{j})-E_{\gamma}]}$ for $j\in\{1,2\}$. Note that we take the positive root such that $p_{j}^{\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,Miller et al. (2003a); Vaníček et al. (2005) the dividing surfaces according to choice A in Table 1 are not close to the turning points of the semiclassical instanton, and under a semiclassical analysis there is no relation between the two choices.

Using Eqs. (6), (20) and (26), we find the semiclassical limit of $C_{\mathrm{ff}}(0)$ to be

[TABLE]

It is important to note that the terms proportional to $K_{\ell}^{2}$ and $K_{r}^{2}$, which appear in $C_{\mathrm{dd}}(0)$, completely cancel out in the semiclassical limit of $C_{\mathrm{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, $\dot{K}_{\gamma}=i\frac{\partial{K_{\gamma}}}{\partial{\tau_{\gamma}}}$ and $\dot{E}_{\gamma}=i\frac{\partial{E_{\gamma}}}{\partial{\tau_{\gamma}}}$, as

[TABLE]

and

[TABLE]

The combination of these two equations and Eq. (13) gives the energy variance of $C_{\mathrm{dd}}$ as

[TABLE]

A similar analysis gives $\Delta H_{\mathrm{ff}}$ in terms of the semiclassical quantities:

[TABLE]

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_{\mathrm{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_{\ell}^{2}$ only cancels to first order in $\hbar$ and cannot be neglected when compared with the much smaller $K_{\ell}K_{r}$ term. This problem is avoided by choosing dividing surfaces according to the standard prescriptions in Table 1, for which $K_{\ell}$ 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_{\ell}=K_{r}$. In the following, we shall use this fact to explain the success of the QI method for symmetric systems.

III.2 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.Richardson (2018b) 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 $\rho$ 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_{\ell}=E_{r}$ and hence $p^{\ell}_{1}=p^{r}_{1}=p_{1}$ and $p^{\ell}_{2}=p^{r}_{2}=p_{2}$. The flux-flux correlation function [Eq. (27)] thus simplifies to

[TABLE]

and

[TABLE]

In order to study the QI approximation based on $\Delta H_{\mathrm{dd}}$, we define dividing surfaces which obey Eq. (10) in the semiclassical limit [see also Eq. (26)], which implies that $K_{\ell}=K_{r}$. The semiclassical limit of the energy variance $\Delta H_{\mathrm{dd}}$ [Eq. (30)] becomes identical to that of $\Delta H_{\mathrm{ff}}$ given in Eq. (33).

Eqs. (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)Gutzwiller (1990)

[TABLE]

and

[TABLE]

Eq. (32) becomes

[TABLE]

UsingRichardson (2016a)

[TABLE]

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

[TABLE]

which matches Miller’s original SCI expression Miller (1975) in the one-dimensional case, bearing in mind that the total action of the instanton periodic orbit is given by $S_{\ell}+S_{r}$. This equation is equal to the ImF instanton expression,Affleck (1981) which can be written in a number of equivalent waysAlthorpe (2011); Richardson (2018b) including the ring-polymer formulation.Richardson and Althorpe (2009) 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 $\beta\hbar/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 Richardson (2018b) 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

[TABLE]

where the projection operators are written in terms of the Heaviside step function, $\theta$, as

[TABLE]

such that

[TABLE]

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, $\braket{x_{2}}{\hat{U}_{\ell}(-i\tau_{\ell})}{x_{1}}$ will be dominated by the left-bouncing minimum-action path and $\braket{x_{2}}{\hat{U}_{r}(-i\tau_{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 either both bounce 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.Richardson (2016a, 2018b) We therefore propose to approximate the rate as

[TABLE]

where

[TABLE]

which is a general (non-symmetric and complex) correlation function whose form permits the imaginary time to be different for the two propagators, although their sum, $\tau_{\ell}+\tau_{r}=\beta\hbar$, 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.Richardson (2016a, 2018b) It is absolutely necessary to project onto left and right paths when using a non-symmetric split of the imaginary time to avoid finding the spurious minimum-action paths of a right-bouncing path in time $\tau_{\ell}$ 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. Richardson (2018b)

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

[TABLE]

where the projected matrix elements are given by

[TABLE]

and their time derivatives calculated using $\dot{\rho}_{\gamma}=i\frac{\partial{\rho_{\gamma}}}{\partial{\tau_{\gamma}}}$ and $\ddot{\rho}_{\gamma}=-\frac{\partial^{2}\rho_{\gamma}}{\partial\tau_{\gamma}^{2}}$. The PQI rate expression is then defined identically to the 2OCE approximation [Eq. (18)] except that $C_{\mathrm{ff}}$ is replaced by $C_{\mathrm{ff}}^{\mathcal{P}}$ throughout.

By following a similar approach as in the previous section we obtain the semiclassical limit of $C^{\mathcal{P}}_{\mathrm{ff}}(0)$ equal to that of Eq. (32), and of $\Delta H_{\mathrm{ff}}^{2}$ 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.2. This result is actually unsurprising, as the derivation of SCI follows similar lines of reasoning in Ref. Richardson, 2018b.

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 $\tau_{\ell}$ (keeping $\tau_{\ell}+\tau_{r}=\beta\hbar$ fixed) until $\dot{C}^{\mathcal{P}}_{\mathrm{ff}}(0)=0$. However, the value of $\tau_{\ell}$ 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.Wolynes (1987) 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 Mills et al. (1997); Richardson and Althorpe (2009); Hele and Althorpe (2013) and most recently in Ref. Thapa, Fang, and Richardson, 2019 for nonadiabatic rates.

The PQI method is also applicable to multidimensional problems and can be efficiently computed using 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^{\prime}$ variable in Eq. (40) and, for example, only paths for which $x^{\prime}<x_{0}$ should contribute to $\rho_{\ell}$. In this work, however, we implement the projection by integrating numerically over the allowed range of $x^{\prime}$ 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 2. 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_{\ell}\gg K_{r}$. The good agreement is also seen in Table 3, 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 formulae 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 analysing the QI results. The PQI method predicts a rate of $4.34\times 10^{-12}$ a.u. without optimization of the imaginary-time split (i.e.  using the semiclassical ratio of $\tau_{r}/\tau_{\ell}\approx 0.258$), but when it is optimized to $\tau_{r}/\tau_{\ell}\approx 0.298$, the rate prediction is $4.85\times 10^{-12}$ a.u., which is even closer to the exact result (c.f. caption Table III).

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, $\alpha$, and for higher temperatures than that considered in Tables 2 and 3.222The error in the QI rates is slightly larger than those reported in Ref. Miller et al., 2003a because they used an ad hoc correction to the energy variance for their calculations, but this does not affect our conclusions. To generate these results, we have generally favoured 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.Ceotto, Yang, and Miller (2005) 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,Richardson (2016b, 2018b) 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;Ceotto, Yang, and Miller (2005) 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 $\tau_{\ell}$ is not optimized but simply chosen to match the imaginary time for the left bounce in the instanton trajectory.

Close to the cross-over temperature, $T_{\mathrm{c}}=\hbar\omega_{\mathrm{b}}/{2\pi k_{\mathrm{B}}}$ (where $\omega_{\mathrm{b}}$ is the absolute value of the imaginary barrier frequency),Richardson (2018b) 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 $C^{\mathcal{P}}_{\mathrm{ff}}(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 Wolynes (1987); Thapa, Fang, and Richardson (2019) 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.Richardson (2018a) RPMD rate theoryCraig and Manolopoulos (2005b) also dominantly samples instanton paths even for asymmetric systemsRichardson and Althorpe (2009) 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.Habershon et al. (2013)

VII Acknowledgements

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 behaviour of PQI for the case of the free particle, which can be treated analytically. The exact propagator for the free particle is well knownFeynman and Hibbs (1965)

[TABLE]

If we substitute $t\rightarrow-i\tau_{\gamma}$ in Eq. (46) and then insert Eq. (46) into Eq. (45), we find

[TABLE]

The integrals over $x^{\prime}$ can be done analytically with standard resultsGradshteyn and Ryzhik (2000) to obtain

[TABLE]

where $\Phi(z)$ is the Gauss error function, and

[TABLE]

Note that, as required from Eq. (41), $\rho=\rho_{\ell}+\rho_{r}$.

The free-particle delta-delta correlation function is

[TABLE]

and the stationarity condition [Eq. (10)] for the dividing surfaces, equivalent to setting $\partial C_{\mathrm{dd}}(0)/\partial x_{\gamma}=0$, implies $x_{1}=x_{2}$. Taking $x_{0}=x_{1}=x_{2}$, we find $\rho_{\ell}=\rho_{r}=\rho/2$.

For $x_{1}=x_{2}$, the exact quantum flux-flux correlation function for the free particle is Miller, Schwartz, and Tromp (1983)

[TABLE]

Quantities required for the QI rate [Eq. (14)] are thus $\Delta H_{\mathrm{dd}}=\sqrt{2}/\beta$ and $C_{\mathrm{ff}}(0)=(\pi\hbar^{2}\beta^{2})^{-1}$. The resulting QI rate,

[TABLE]

is a factor of $\sqrt{\pi/2}$ from the correct free-particle result, $kQ_{\mathrm{r}}=(2\pi\hbar\beta)^{-1}$, leading to an overestimation of the rate by 25$\%$.Miller et al. (2003b)

For the 2OCE rate (18), the relevant quantity is instead $\Delta H_{\mathrm{ff}}=\sqrt{6}/\beta$, giving the rate

[TABLE]

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

Likewise, we can calculate the projected correlation function for the free particle using Eqs. (44) and (48) with the substitution $\tau_{\gamma}\rightarrow\tau_{\gamma}+it$ to obtain

[TABLE]

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

[TABLE]

and

[TABLE]

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

[TABLE]

which underestimates the exact rate by about 37%.

It is possible to directly integrate $C^{\mathcal{P}}_{\mathrm{ff}}(t)$ over $t$ to yield the exact rate for this system. Then using Eq. (42), we get $kQ_{\mathrm{r}}=(2\pi\hbar\beta)^{-1}$, which is the exact rate for this system. Miller et al. (2003a) 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 Miller et al. (2003a); Ceotto, Yang, and Miller (2005); Hansen and Andersen (1996) 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. EckartEckart (1930) 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

[TABLE]

where comparison with Eq. (19) gives $\xi=-e^{2\pi x/l}$, $A=V_{0}(1-\alpha)$, $B=V_{0}(1+\sqrt{\alpha})^{2}$ and $l=\pi/a$. Following Eckart, we define

[TABLE]

and then rearrange into the form

[TABLE]

Note that if $E<0$, the root in Eq. (59a) is taken such that $\operatorname{Im}\bar{\alpha}<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. (60). Klein Klein (1933) showed that solutions to this differential equation can be obtained using the ansatz

[TABLE]

with

[TABLE]

However, $F[a_{0},b_{0},c_{0},\xi]$ is only one possible solution to the hypergeometric differential equation, and instead we choose to use Forsyth’s solutions Forsyth (1948) called III and XIII to give

[TABLE]

The former is only a solution if $E>A$ and the second if $E>0$. In order to find the $\xi\rightarrow-\infty$ limit, we first apply Abramowitz and Stegun’s formula 15.3.7,Abramowitz and Stegun (1972) and after premultiplying by a constant we find the desired asymptotic limits for $u_{1}$ to be

[TABLE]

where

[TABLE]

and $k_{\bar{\alpha}}=2\pi\bar{\alpha}/l$ and $k_{\bar{\beta}}=2\pi\bar{\beta}/l$. The asymptotic limits of the second wavefunction are (using Abramowitz and Stegun 15.3.9)Abramowitz and Stegun (1972)

[TABLE]

where

[TABLE]

Note that Eckart’s results for reflection and transmission are given by $|B_{2}/A_{2}|^{2}$ and $|\bar{\beta}/\bar{\alpha}|/|A_{2}|^{2}$ given a particle incident on the left.

Finally, we normalize the wavefunctions

[TABLE]

such that

[TABLE]

Using

[TABLE]

we find the normalization constants should be chosen as

[TABLE]

where $k_{\bar{\alpha}}^{\prime}=\frac{\partial{k_{\bar{\alpha}}}}{\partial{E}}=k_{\bar{\alpha}}/(2E)$ and $k_{\bar{\beta}}^{\prime}=\frac{\partial{k_{\bar{\beta}}}}{\partial{E}}=k_{\bar{\beta}}/(2(E-A))$. 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 $B_{2}/k_{\bar{\alpha}}^{\prime}+C_{1}^{*}/k_{\bar{\beta}}^{\prime}=0$.

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

[TABLE]

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.