Instanton theory for Fermi's golden rule and beyond.

Instanton theory provides a semiclassical approximation for computing quantum tunnelling effects in complex molecular systems. It is typically applied to proton-transfer reactions for which the Born-Oppenheimer approximation is valid. However, many processes in physics, chemistry and biology, such as electron transfers, are non-adiabatic and are correctly described instead using Fermi's golden rule. In this work, we discuss how instanton theory can be generalized to treat these reactions in the golden-rule limit. We then extend the theory to treat fourth-order processes such as bridge-mediated electron transfer and apply the method to simulate an electron moving through a model system of three coupled quantum dots. By comparison with benchmark quantum calculations, we demonstrate that the instanton results are much more reliable than alternative approximations based on superexchange-mediated effective coupling or a classical sequential mechanism. This article is part of the theme issue 'Chemistry without the Born-Oppenheimer approximation'.

Introduction

The Born–Oppenheimer approximation is widely used to simplify the Schrödinger equation for molecular systems by expressing the equations of motion in terms of a single adiabatic potential-energy surface (PES). However, there are many molecular processes for which this adiabatic approximation is not valid, typically due to two PESs becoming close in energy. In these situations, alternative approaches are necessary.

One way of going beyond the Born–Oppenheimer approximation is to represent the total Hamiltonian in the basis of diabatic electronic states [1]. If the coupling between the diabatic states is small, perturbation theory can be used to obtain solutions correct to a given order. Typically, the leading-order approximation is second order in the coupling, and the resulting expression is commonly called Fermi’s Golden Rule (FGR) [2]. Within this theory, the rate for a transition from an initial state to a final state is [3-5] which is proportional to the square of the coupling between the states, , and by virtue of the delta function, enforces energy conservation. This expression can alternatively be written in terms of quantum-mechanical propagators or Green’s function operators [6,7] and is also known as ‘the first Born approximation’ in scattering theory [8].

FGR is very general and can in principle be applied to study any physical system that evolves between a pair of weakly coupled initial and final states. It has found applications in various fields, including spectroscopy [9,10], scattering processes such as X-ray diffraction [6,8], scanning tunnelling microscopy [11], and the kinetics of nuclear [2] and electron-transfer [12,13] reactions.

Marcus theory can be recovered from equation (1.1) by taking into account the initial thermal distribution of particles and replacing the quantum states with a classical phase-space average under the assumption of parabolic free-energy surfaces [14]. An important consequence of this theory was the prediction and later experimental verification of the inverted regime [15]. However, these classical theories neglect quantum tunnelling, and many quantum-mechanical methods developed to extend Marcus theory rely on the assumption of harmonic PESs, for which the analytic solutions to the Schrödinger equation are known [12]. These approaches thus constitute uncontrolled approximations which may not always reflect the reality of complex anharmonic molecular systems. Instead, we require a formalism which provides a rigorous asymptotic approximation to the FGR rate without the necessity of first solving the Schrödinger equation, even for the reference Hamiltonians.

Semiclassical instanton (SCI) theory [16,17] provides a framework for describing multidimensional tunnelling effects in complex molecular systems. It is based on the path-integral formulation of quantum mechanics [6] but avoids summing over all possible paths by identifying the optimal tunnelling pathway, which is called the ‘instanton’. The instanton pathway is defined as the path which makes the action stationary. It is represented by a periodic classical orbit in imaginary time, or equivalently one which obeys Newton’s equations of motion on an upside-down PES [16]. The use of imaginary time to simulate tunnelling processes can be justified by considering a particle crossing a potential-energy barrier, , with energy lower than the barrier top. In order to pass through the region where , the kinetic energy must take negative values, which implies that the particle will have imaginary momentum and can therefore be said to travel in imaginary time.

The ring-polymer formulation of instanton theory provides computationally efficient algorithms for locating the optimal tunnelling pathways on a single Born–Oppenheimer PES and evaluating the rate [17-22]. Because the method is based on a unique optimal pathway and does not rely on sampling a statistically large number of path configurations, it can be applied to complex molecular systems which require expensive electronic-structure calculations.[1] In addition to calculating rates accurately, the instanton pathway also offers clear mechanistic insight into the tunnelling process. For instance, quantitative information on the contribution from each atom can be extracted such that kinetic isotope effects can be easily predicted and understood [26-28]. It has, for instance, been applied to gas-phase reactions [29-31], tunnelling in molecules and clusters [28,32-35], enzymatic reactions [36], surface diffusion processes [37,38] and water formation reactions [39].

An extension of the original instanton theory has been derived which provides a semiclassical approximation to FGR rates of non-adiabatic reactions and predicts very good results for benchmark model systems [40-42]. Unlike the Zhu–Nakamura formulae [43] (which extend Landau–Zener theory [44,45] to include tunnelling along a one-dimensional reaction coordinate), instantons constitute a multidimensional tunnelling theory. They are thus expected to provide more accurate results than the one-dimensional approaches (assuming the reaction takes place in the golden-rule limit), although at an enhanced computational cost as information is required along the instanton path and not just at a single point. The ring-polymer discretization of the instanton trajectory again provides a computationally efficient algorithm in both the normal [46,47] and inverted regimes [42], and we present the working equations in §3.

For many chemical processes, however, it is necessary to go beyond the golden-rule approximation. This is particularly clear for systems where the initial and final states are both coupled to a third ‘bridge’ state but not to each other. This process is sometimes referred to as the superexchange mechanism [48-50] or the through-bond coupling mechanism [51]. Examples include bridge-mediated electron transfer in bacterial photosynthetic reaction centres [52], hole migration in DNA [50], CO-haemoglobin recombination [12] and two-photon processes such as Raman spectroscopy [10]. Such processes are described by higher-order perturbation-theory, leading to the fourth-order rate constant [5] This describes a concerted mechanism from to via intermediate bridge states [49].

Just as for the second-order FGR case, equation (1.2) can only be evaluated if one first solves the time-independent Schrödinger equation for the uncoupled problem in order to obtain the energy levels , and . For complex molecular systems, this is practically impossible without making further approximations. Even if reasonable approximations are obtained for the energy levels, the calculation of the sum over the intermediate states can be challenging. This is because the individual terms in the sum, , may be either positive or negative, and a large number of the terms are similar in magnitude but differ in the sign. This makes it difficult both to evaluate the sum numerically and to extract a simple mechanistic interpretation from the expression.

In this paper, we derive a new fourth-order instanton approximation as an extension of the derivation in the golden-rule limit. This provides reasonably accurate numerical results as well as a simple mechanistic interpretation for the reaction. We start by reviewing the existing golden-rule framework before tackling the fourth-order case in §4.

Golden-rule instanton theory

In this section, we outline the derivation of golden-rule instanton theory starting from the standard expression for Fermi’s golden rule.

Fermi’s golden rule

We consider the total Hamiltonian where and are the diabatic electronic states of the reactant and product.

The uncoupled Hamiltonians are typically given in terms of the PESs, , by where . This describes a system with nuclear degrees of freedom with coordinates and conjugate momenta which have been mass-weighted such that the associated mass of each degree of freedom is . The electronic states are coupled by , which is assumed to be small and slowly varying.

In the standard perturbation theory approach, we assume that we can solve the time-independent Schrödinger equation for the unperturbed Hamiltonians Each set of vibrational states (with quantum numbers or ) forms a complete basis of the nuclear Hilbert space such that . The vibronic states and form a basis for the total Hamiltonian and are coupled by .

From equation (1.1), taking a Boltzmann average over the initial vibronic states , for which the reactant partition function is , and summing over all final states , the second-order thermal rate at inverse temperature is given by , where To derive equation (2.4b), we first replaced by . This is valid for any value of due to the presence of the delta function, which forces the two energies to be equal. The reason for this manipulation will become apparent later. Next, we replaced the delta function by its Fourier integral representation. This allows us to rewrite the rate expression in a basis-independent form without reference to the eigenstates of the uncoupled systems and paves the way for a path-integral formulation. The quantum-mechanical rate can thus be defined by with the correlation function This is related to the general flux correlation function formulation of rate theory [53-55] (although subtly different due to the asymmetric shift introduced by [18,24]) and is the starting point for a number of path-integral simulations of the golden-rule rate [25,41,46,56-60]. Note that can be thought of as an analytic function of and that modifying is thus equivalent to deforming the integration contour in the complex plane. Thus, although explicitly depends on , the Cauchy integral theorem [61] ensures that is independent of , which will allow us in the following sections to pick the most convenient value.

Semiclassical instanton theory

We have shown how the FGR rate can be expressed in terms of the second-order correlation function, . In this subsection, we derive an SCI approximation to the correlation function and hence to the FGR rate itself.

First, we write the correlation function in the basis of position states where the quantum-mechanical propagators, , give the probability amplitude for a particle starting at to travel to in complex time (the real part of corresponds to imaginary time and the imaginary part to real time). We will see later that the time integral over the correlation function is dominated by its value at , and thus the propagators can be approximated using the imaginary-time generalization [62] of the van-Vleck semiclassical propagator [63] Here, the action is defined as for the imaginary-time classical trajectory which solves , i.e. Newton’s equation of motion on the upside-down PES [62], with the boundary conditions and . The prefactor is given by the determinant of the second-derivative matrix of the action:

The remaining integrals in equations (2.5) and (2.7) are evaluated using the steepest-descent approximation, which, for the general case given below, is derived by expanding the exponent to second order around its stationary point and the prefactor to lowest order, where is a vector of dimension , and are arbitrary analytic functions, , is the solution of and we have assumed that is positive definite and that . For our purposes, the integration range is typically over the whole real axis, but the expression is valid as long as is inside the range. In the following, we shall omit the tilde; whenever appears outside an integral, we imply that it is evaluated at the stationary point . It can be proved that the steepest-descent approximation is an asymptotic expansion of the integral in and is thus exact in the limit [64]. The van-Vleck propagator itself can be derived as a steepest-descent approximation to the exact path-integral propagator [63,65].

We then arrive at the SCI approximation by replacing the propagators with their semiclassical limits (equation (2.8)) and performing the integrals over space and time simultaneously with the method of steepest descent (equation (2.11)). Typically, a stationary point can be found with real-valued , and , which corresponds to . The final expression for the rate is thus where and all expressions (including ) are evaluated with , and at the values that make the total action stationary. The prefactor includes the determinant and the negative sign appears as a consequence of the Cauchy–Riemann equations, . As we show in appendix A, an equivalent derivation also holds in the inverted regime. For consistency, we choose to evaluate the reactant partition function, , using an equivalent semiclassical limit, which typically corresponds to the standard quantum-mechanical harmonic approximation [18].

The stationary-point conditions, and , give important information about the instanton. In particular, because (and equivalently with double primes), momentum is conserved at the hopping points. Also, if we identify and , we find that because , energy must be conserved. The hopping points must therefore be located on the crossing seam, where , but not necessarily at the minimum-energy crossing point (MECP).

In previous work, we have applied the instanton approach to various model systems including asymmetric system-bath models of electron transfer [47] and photodissociation in anharmonic potentials [42]. The results compared favourably with quantum-mechanical benchmarks and were seen to be much more accurate than alternative methods such as classical Marcus theory and the cumulant expansion [66].

Typically, for molecular systems in the gas phase and even many systems in the condensed phase, we have to take into account that in addition to the vibrational degrees of freedom the molecule also possesses external degrees of freedom comprising translations and/or rotations. The presence of these external degrees of freedom causes some of the eigenvalues of to be zero, preventing a straightforward application of the steepest descent formula. We can, however, project them out and account for these degrees of freedom by the standard expressions for the partition functions for global translation and rotation.

In order to account for the external degrees of freedom we first consider the free-particle action given by . Because , it can readily be seen that translations and rotations manifest themselves in in the form of eigenvalues equal to . A similar argument shows that has zero eigenvalues and eigenvalues equal to resulting from the external degrees of freedom. The contributions of the translations and rotations can thus be separated by removing these eigenvalues in the calculation of the determinants and . Finally, under the assumption that rotations and translations do not couple to the vibrations, one multiplies the rate equation (2.12) by the translational and rotational partition function corresponding to the MECP or averaged along the instanton [18].

Application to the separable case

If we apply the instanton formulation of the previous subsection to a separable Hamiltonian, we obtain a closed-form formula for the golden-rule rate. In [41], a fully classical rate formula was presented based on the knowledge of the minimum-energy crossing point. We go a step further and employ a separable model which includes tunnelling along an arbitrary one-dimensional reaction coordinate and quantized vibrations in the perpendicular modes in an analogous way to Eyring’s transition-state theory (TST) [67]. The formula we obtain is similar but not identical to those commonly used in the literature [7,12,68-70]. The key approximation here neglects frictional effects (i.e. coupling between the reaction coordinate and the other modes) which would otherwise cause the instanton path to curve and would require a full-dimensional numerical solution as described in §3.

The most important configuration for reactions in the golden-rule limit is the minimum-energy crossing point, , defined by the minimum of subject to the constraint , which generalizes the transition-state concept to non-adiabatic rate theory. At this point, the potential is and the two gradients, and , are either parallel or antiparallel, corresponding to the inverted or normal regime. The Hessians are ; note that in general it will not be possible to find a normal-mode transform which simultaneously diagonalizes and .

The potential-energy surface is assumed to be separable around the MECP and we use a coordinate system consisting of a reaction coordinate, , which points in the direction of the MECP gradient, , and a set of coupled harmonic vibrational modes, , for the degrees of freedom perpendicular to the reaction coordinate. The origin of this coordinate system is taken to be the MECP itself. In many molecular systems, there are also modes with zero frequencies corresponding to translations and rotations, which will be designated by . The most common case is a nonlinear molecule of atoms in vacuum, which has degrees of freedom in total, of which are zero modes. Three of the coordinates thus determine the location of the molecular centre of mass, and three of them are Euler angles mass-weighted according to the moments of inertia. Here, we will neglect Coriolis terms [71] and keep the moments of inertia constant using their values at the MECP.

The diabatic potentials are thus approximated by where and is a matrix which projects out the reaction coordinate and zero modes corresponding to translations and rotations [70,72]. In general, the two curves can take an arbitrary form as long as they cross. In the simplest case, the potential along the reaction coordinate is assumed to be linear: , where is the slope of the PES at the MECP. One level higher is the harmonic approximation, , where and .

Given the form of equation (2.14), the action is separable, and the parts corresponding to the zero modes and harmonic vibrational modes can be written as analytic functions [6]: where , and the matrices are defined by . We have employed the orthogonal coordinate transform and and likewise for coordinates (and later also for ).

The action along the reaction coordinate, is defined analogously to equation (2.9). In the case of a linear approximation, , and for the harmonic approximation, . However, in general, the action can be evaluated using a Legendre transform where the line integral in equation (2.16c) is along the path from one end point to the turning point and back to the other end point. The sign function, , is necessary for treating the inverted regime where [42]. It is thus possible to evaluate all quantities by one-dimensional numerical integration. The total action along the reaction coordinate only is , which is the obvious generalization of the total action of §2(b).

The stationary point is clearly located where and . The remaining variables , and are determined by considering the one-dimensional problem along the reaction coordinate only. At the stationary point, we must have as well as and likewise . These equations enforce energy and momentum conservation and can be solved only at the crossing point, implying that .

The prefactor terms are given by [18,63]

Because in many cases mixed derivatives of position and time vanish at the stationary point, the second-derivative matrix separates into blocks such that the determinant , where has the form of equation (2.13) for the action and In general, and have different normal modes and cannot therefore be simultaneously diagonalized. It is therefore not possible to simplify the expression further. Nonetheless, they can at this point be easily evaluated given the Hessians of each state.

The contribution from the zero modes is , which is equivalent to the standard prefactor of a classical configurational integral. Thus, although integration over the zero modes cannot be performed by steepest descent, it can be substituted by the classical partition functions instead. This gives a rate (within the separable approximation) where the transition-state partition function is . Here, the translational and rotational partition functions are given by the usual expressions and are evaluated at the MECP. The reactant partition function, , takes the standard form within the (classical) rigid-rotor and (quantum) harmonic-oscillator approximations and takes account of whether a bimolecular or unimolecular reaction is studied. The vibrational partition function of the MECP is given by and will be discussed later in this section.

First, we study the effect of tunnelling on the rate within the linear approximation, for which the stationary point is found at and leading to where we have defined and the prefactors are and .

In the classical limit, we obtain the non-tunnelling harmonic transition-state theory (hTST) rate [12,68] which is similar in spirit to Eyring’s TST but applicable in the golden-rule limit.

By including the tunnelling effect, the rate is effectively multiplied by a Gamow-like factor:

Although we derived this formula using the instanton method, it also happens to be the exact quantum solution for the linear model [41]. The linear model predicts tunnelling at an energy of . This approach clearly breaks down if as this would correspond to an unphysical initial state. In this case, a more appropriate model is definitely required for .

It is, however, not much more complicated to employ the separable harmonic approximation, which leads to [12,41] where and . It has a maximum at . From this, we obtain an instanton rate of The tunnelling factor thus differs considerably from that used in the linear approximation. We show in figure 1 that using the linear approximation can lead to predictions for the tunnelling effect which are orders of magnitude too large if the system deviates from the linear form. The harmonic model tested uses reasonable values for the gradients of the slopes in the inverted regime on the order of those reported in [70]. Here, we have assumed the same frequency for the reactant and product states, but it could be further generalized to treat asymmetric frequencies at the cost of using more complicated formulae [47,73,74].

Figure 1

Golden-rule rates obtained for the separable harmonic system using the linear (equation (2.24)) or instanton (SCI) approximation (equation (2.26)) are compared with quantum-mechanical results. Two versions of the WC approximation (equation (2.27)) are also presented which differ in the value of used as indicated. The fixed parameters are , at 300 K. The quantum-mechanical golden-rule result (equation (2.5)) is defined by integrating the correlation function until the first plateau. The inset shows the percentage error made by the instanton approximation over the same range of . (Online version in colour.)

We also compare with another commonly used formula, the so-called weak-coupling (WC) approximation, which is defined by [70,72,75]

where the transmission probability is given in terms of the Airy function as [75,76] and was defined in equation (2.22). Two suggestions have appeared in the literature [77], either to choose as the minimum potential energy of the reactants, or the zero-point energy of the reactants. Setting exactly recovers the expression in equation (2.24). The WC approximation is valid when the integrand is dominated by energies near the crossing point, but there is no formal justification for the approach in the deep tunnelling regime, regardless of how is defined, hence the deviation from exact results observed in figure 1.

We now return to the discussion of the vibrational partition function. The formula for in equation (2.20) differs from that typically used in the literature. Probably, the most common approach is to quantize the vibrations of an effective Hessian defined by [70,72].

The two methods are in agreement in the classical limit as we can prove by expanding the trigonometric functions to first order in their arguments to obtain where . If we take the solution of from the linear approximation (which is also valid at high temperature for an arbitrary system) [41], we recover .

However, at low temperatures, the commonly used approximation [70], can only be derived from equation (2.20) in the limiting case that . In conclusion, the usual formula (equation (2.29)) is only rigorously justified in the classical limit or in the case that the difference between the two Hessians is negligible. However, for consistency with the separable assumption made in equation (2.14), we recommend that one should use the value of determined from the stationary point and the general formula, equation (2.20).

Ring-polymer formulation

In this section, we present the working equations for the ring-polymer formulation of the golden-rule instanton. In particular, we go beyond previous work [42,46,47] in three aspects: we present an optimization algorithm using only half of the ring polymer; we obtain the rate directly from the eigenvalues of the ring-polymer Hessian; and we explicitly treat translational and rotational degrees of freedom. This makes the algorithm as similar as possible to the standard ring-polymer instanton method within the Born–Oppenheimer approximation [18].

Ring-polymer rate

As described in [42,46], we discretize the instanton path and its associated action in the form of a ring polymer. First, the instanton is located using a saddle-point search in the extended ring-polymer space as described in §3(b). Then, instead of evaluating the rate via equation (2.12), which is defined in terms of continuous trajectories and van-Vleck propagators, it is more natural to obtain the result directly from the ring-polymer picture. This approach is in direct analogy to the standard Born–Oppenheimer ring-polymer instanton method [17-20,78,79].

The ring-polymer rate expression can be readily derived by applying the Trotter formula to split each of the propagators in equation (2.7) into intervals of imaginary time, , and inserting a complete set of position states between each interval before applying the steepest-descent approximation. We thus discretize the instanton orbit into ‘beads’ which constitute replicas of the system along the pathway. In order to draw a direct connection to the application of the method to molecular systems, we introduce as the three-dimensional atomic coordinate vector for atom of bead with corresponding mass . The atomic vectors for bead can be combined into the -dimensional coordinate vector , where and is the number of atoms. This results in the following expression for the correlation function at :[2] where is the combined vector of all bead positions, and the action is where cyclic indices are implied such that , forming a closed ring. Equation (3.2) is defined such that in the limit the value of the ring-polymer action at its stationary point is equal to the SCI action, . Equivalently to the derivation in §2(b), we seek to obtain the instanton rate expression by simultaneous steepest-descent integration about the stationary point of in and . However, before this can be done, we must account for the zero-frequency modes corresponding to global translations and rotations of the ring polymer. Note that, in contrast to the standard Born–Oppenheimer instanton, there is no zero mode corresponding to permutations of the beads here.

In analogy to [18], we will write the ring-polymer instanton rate as a product of translational, rotational and vibrational contributions where , which is equivalent to in the normal but not in the inverted regime (where ), and is the ring-polymer potential of the instanton configuration. The number of beads on the reactant and product states can be chosen arbitrarily and independently of each other, which gives rise to an unappealing feature. Namely, even though the total rate is independent of the ratio as long as the instanton is converged with respect and , the values of the rotational and vibrational partition functions can change depending on the user’s choice for the bead split when the standard definitions for the ring-polymer partition functions are employed. We can resolve this issue by adopting a bead-dependent mass-weighting scheme that renormalizes the individual partition functions and renders them independent of the chosen bead distribution.

We define the bead-dependent masses This mass-weighting scheme does not alter the standard ring-polymer expression for , as the total mass of the ring polymer remains equal to , where is the total mass of the molecule. The moment-of-inertia tensor for the ring polymer, however, is affected by the mass-weighting where is the identity matrix, is the outer product and we have assumed that the centre of mass has been placed at the origin, i.e. . Note that in equation (3.5) the mass is not only different for different atoms but also for different beads, in line with our definition in equations (3.4). Given , the rotational partition function is then defined by the standard equations for linear or nonlinear configurations (e.g. eqn (60) in [18]).

The vibrational partition function is defined by where are the eigenvalues of the mass-weighted ring-polymer Hessian defined by , with the gradient with respect to bead positions and imaginary time combined ] and the mass matrix containing the bead-dependent masses defined in equations (3.4). In analogy to the coordinate vector of bead , we defined the corresponding mass vectors containing the masses associated with all degrees of freedom. Here, we also choose to weight the imaginary-time coordinate by to ensure all entries of the ring-polymer Hessian have the same unit. The zero eigenvalues corresponding to global translations and rotations are excluded from the product in equation (3.6). The eigenvalues associated with the single imaginary mode in the normal regime, or the imaginary modes occurring in the inverted regime (see appendix A and [42]), enter the product with their absolute magnitude. Finally, there is a factor , which accounts for the masses and imaginary-time intervals in the prefactor from equation (3.1), as well as the mass-weighting. It is given by Note that this factor is equal to 1 for the mass-weighting scheme introduced above if the bead ratios are chosen such that . We emphasize that different choices of mass-weighting schemes if implemented correctly (i.e. taking the factor into account) do not affect the total rate, but only the relative values of the vibrational, rotational and translational contributions. For a molecule with no translational or rotational degrees of freedom, the ring-polymer instanton rate (equation (3.3)) reduces to equation (2.12) as the number of beads tends to infinity.

Optimization algorithms

In this subsection, we will discuss various numerically efficient algorithms that can be used to locate golden-rule instantons within the ring-polymer formalism. As has been mentioned earlier, this is equivalent to finding a saddle point of the ring-polymer action (equation (3.2)). In the normal regime, as for the Born–Oppenheimer instanton, this saddle point is of index one, i.e. the second-derivative matrix of the action exhibits one negative eigenvalue, and thus the standard transition-state optimization routines can be applied, typically quasi-Newton algorithms [80-82]. The key difference is that we must also simultaneously optimize along with .

Analogously to the Born–Oppenheimer instanton [18,20,83], the trajectories of the golden-rule instanton typically fold back on themselves [46]. In these cases, we can approximately halve the computational effort required for evaluating the action and its derivatives. We may therefore locate the instanton by searching for a single-index saddle point of the half ring-polymer action, which is defined such that for the instanton configuration and has been redefined for the half ring polymer. Note that the bead labelling has changed compared with equation (3.2), since cyclic indexing does not apply to the half ring polymer and there is only a single hopping bead (figure 2).

Figure 2

Illustration of how a half ring polymer consisting of six beads (dashed box) relates to the corresponding 10-bead full ring polymer. Here, and . Note that the configuration of beads shown here is typical for an instanton in the inverted regime, although the numbering scheme is identical in the normal regime. (Online version in colour.)

When optimizing an instanton for a new system, one can use the fact that at high temperatures the instanton will be compact and located near the MECP, which can itself be found with standard methods [84-90]. A good initial guess for the value of is given by the classical, short-time limit of the action (see discussion above equation (2.21)) [41]. It is therefore relatively easy to optimize an instanton at high temperature, which can then serve as a starting point to slowly cool down the system. This process is accompanied by the gradual expansion of the instanton.

In contrast to the normal regime, in the inverted regime of an exoergic reaction, the ring-polymer instanton corresponds to a saddle point of index .[3] In the inverted regime, we can also employ the half ring-polymer approach (equation (3.8)), though here the index of the required saddle point changes to . From figure 2, it can be seen how the half ring polymer emerges from the full ring in a configuration typically assumed in the inverted regime.

Due to the fact that in the inverted regime one is no longer seeking a single-index saddle point, different optimization algorithms are required. Conveniently, we can make use of algorithms originally developed for locating higher-index saddle points in Lennard–Jones clusters [91-93]. These algorithms are based on eigenvector following, where uphill steps are taken in the direction of the eigenvectors corresponding to the lowest eigenvalues and downhill steps are taken in the remaining degrees of freedom. Hence, even if the initial guess does not have the required number of negative eigenvalues, the eigenvector-following algorithm will guide us towards the stationary point of the correct index. In principle, it is also possible to employ standard Newton–Raphson optimizers. These, however, require an initial guess with the precisely right index, since the algorithm cannot change the signs of eigenvalues.

As in conventional geometry optimization, it is often helpful to project out global rotations and translations. This can be done with standard methods [93], treating the ring polymer as a ‘supermolecule’ [18]. It is however necessary to extend the projection matrices into the combined space by appending zeros for the extra degree of freedom.

Application to molecular dissociation

Consider the rate of unimolecular dissociation of a dimer AB defined by the intramolecular pair potentials where is the bond length. We chose physically reasonable parameters similar to those of [94] (although the process studied in that paper is subtly different): , , , , and a reduced mass of . The dimer is initialized in a thermal equilibrium of the bound metastable state and decays via electron transfer into the state, from where it can dissociate.

Quantum-mechanical results were obtained from Fermi’s golden rule by numerically solving the effective one-dimensional Schrödinger equations using quadrature on a polynomial basis for the bound potential and Numerov’s method for the unbound state. In each case, a centrifugal term was added to the potential, where [95,96], to give a rate thermally averaged over vibrational modes for the given angular momentum quantum number and associated partition function . Finally, rates were summed over the angular momentum quantum number, , and weighted by the degeneracy, to give the total thermal rate , where . Figure 3a shows the contributions from the various terms evaluated using the exact quantum expressions, as well as via SCI theory and harmonic transition-state theory. The instanton approximation leads to a slight overprediction of the rates, whereas the hTST results are orders of magnitude too small. In both cases, however, it can be seen that for this diatomic there is little deviation between the predictions based on summing over angular momentum quantum numbers, and those based on the rigid-rotor approximation. In polyatomic molecules, the situation is expected to be even more favourable, as the moments of inertia are typically larger and do not change considerably along the tunnelling pathway. Note that the effective rotational constant in the instanton picture is , whereas for hTST it is given by the value at the MECP, . As the predicted total thermal rate is inversely proportional to this value, the agreement with the quantum results has been improved by using instead of .

Figure 3

Results for the rate of dissociation of the model dimer system defined in equations (3.9). (a) Contributions, , to the thermal rate at from different values of angular momentum quantum number, , where . The dots are evaluated using the semiclassical instanton (SCI) approximation based on an effective one-dimensional system including the centrifugal potential, whereas the lines of the same colour correspond to a rigid-rotor approximation. Results are normalized by the total quantum-mechanical rate . The inset shows the same plot on a logarithmic scale and in addition includes harmonic transition-state theory (hTST) results. (b) Thermal reaction rates computed with SCI, hTST and quantum mechanics. (Online version in colour.)

In figure 3b, results are presented for the golden-rule rate constant obtained by instanton theory and harmonic TST, which can be compared with a benchmark quantum calculation. As expected, there is a plateau in the Arrhenius plot, which is also observed as a manifestation of low-temperature tunnelling effects within the Born–Oppenheimer paradigm [97]. This can be explained within the quantum-mechanical picture, as only the ground state is occupied at low temperatures. In instanton theory, the tunnelling energy tends to zero (using the minimum of as zero on the energy scale) as the temperature is lowered, and tends to a constant, thus explaining why the rate becomes independent of . The semiclassical approximations behind instanton theory clearly remain reasonably accurate even in this extreme regime, with the rates overestimated by only 20%. Note that, even though the instanton energy tends towards the bottom of the reactant well, the theory does not neglect the effect of ZPE, which is instead dealt with in the prefactor terms. In fact, the same thing happens in the Born–Oppenheimer case, for which it is possible to prove that the low-temperature limit of instanton theory is asymptotic to the WKB rate of tunnelling out of the lowest vibrational level [97]. Because of this the instanton results are in good agreement with the quantum calculations in the deep-tunnelling regime and correctly recover the classical limit at high temperatures.

Fourth-order instanton theory

We consider a system that consists of three potentials corresponding to a reactant , bridge and product , arranged as shown in figure 4a. In particular, we exclusively consider the case for which the reaction is in the normal regime and the minimum of the bridge is higher in energy than both the reactant and product minima, but at the point where the reactant and product cross (with energy ), the bridge PES is lowest. We assume that and are coupled to but that is not directly coupled to . This may be rigorously a consequence of selection rules based on symmetry, or only an approximation in which it is valid to neglect the direct mechanism due to the higher-energy crossing point. It is clear that for this problem, we require a fourth-order rate expression, such as equation (1.2), which is given by Dirac in his classic textbook [5]. In this section, we will derive an instanton approximation for this rate.

Figure 4

(a) A simple model system composed of reactant, bridge and product potential-energy surfaces. The vibrational energy levels of each PES are given by thin horizontal lines. The instanton is shown by a thick horizontal line at the optimal tunnelling energy. (b) A schematic illustrating the instanton for an arbitrary set of PESs. The trajectories are shown with bold lines and the PESs with dashed lines, in each case with the colour corresponding to the diabatic state. The label below a trajectory denotes the imaginary time for travelling along it. The hopping points , , and are shown as black dots and labelled above the corresponding point. (Online version in colour.)

Derivation of the fourth-order instanton rate

We write the total Hamiltonian as a sum of an uncoupled system and a small perturbation with where and are defined as before and is the Hamiltonian corresponding to the bridge PES, . The derivation that follows uses a similar approach to that of the second-order rate from §2(a). In addition to equation (2.3), we now also have the eigenvalue equation

The system we are studying evolves from an initial vibrational state to a final vibrational state via the intermediate vibrational states . If we average equation (1.2) over a thermal Boltzmann distribution of the reactant states (with partition function as before), the rate expression for this system becomes Note that we have chosen to write the expression in this symmetric form using the presence of the delta function to freely interchange the energies of the initial, , and final states, .

In order to derive an instanton approximation, we write the delta function and the terms as integrals over real and imaginary time, respectively and Note that equation (4.5b) is only valid if . For our problem, this means that the bridge states have to be sufficiently high above the minimum of both the reactant and product PESs, such that the Boltzmann factor makes the terms with negative values of and negligible. The following expressions are therefore asymptotic to the exact rate in the limit of high bridge energies.

The physical interpretation of this assumption is that the reaction proceeds via a concerted single-step mechanism from to via virtual states, rather than a sequential two-step process that goes from to followed by to [49]. In the context of instanton theory, this implies that the instanton should be located at an energy below the minimum of the bridge, i.e. this is a deep-tunnelling process.

Collecting the energy terms, making use of the eigenvalue equations (2.3) and (4.3) and the additional resolution of identity , we find where the three-time correlation function is This correlation function is an extension of the one from the second-order case (equation (2.6)). In between the propagators for and , we now have to propagate on the bridge, .

Using the same arguments as those for equation (2.5) in §2(a), the overall rate calculated using equation (4.6) is independent of the choice of . As in §2(b), we can write this correlation function in terms of imaginary-time propagators

Next, we invoke the semiclassical approximation and use kernels that correspond to classical trajectories in imaginary time from equation (2.8). When we perform steepest-descent integration over all variables, the four trajectories combine into a periodic orbit, which we identify as the generalized fourth-order instanton with action with each component defined according to equation (2.9). The resulting trajectory is shown in figure 4b. In complete analogy with the golden-rule case, the stationary-point conditions ensure that we have momentum and energy conservation between each segment of the path, leading to a periodic orbit of imaginary time .

The instanton approximation to the rate is therefore where and are defined equivalently to equation (2.10), , where , and similarly for and . The determinant of the second derivatives of the total action is with now denoting the index-one saddle point of . We have again used Cauchy–Riemann identities to write derivatives of in terms of derivatives of . As before, we evaluate the reactant partition function within the semiclassical approximation [18].

In this way, we have considerably simplified the exact quantum-mechanical expression for the three-state case (equation (4.4)), for which carrying out the double sum over the bridge states is tedious and numerically challenging due to the sum containing many terms of similar magnitude and opposite sign. By contrast, instanton theory defines a unique semiclassical pathway which provides a clear picture of the reaction mechanism.

Because the instanton trajectory travels below the bridge, it implies that the bridge states should be considered as virtual states. This is in conflict with the standard assumption that, because the bridge drops below , it should be considered as an intermediate [12] which accepts population transfer in a sequential mechanism [49]. However, as we shall demonstrate below, the instanton formalism predicts the correct behaviour in the low-temperature limit, whereas an expression based on a sequential mechanism does not.

The instanton approach is, however, only applicable if a saddle point of the action exists, which limits the theory to the deep-tunnelling regime. This will always be the situation at low enough temperature, but above a certain crossover temperature, the saddle point will disappear and it is clear that a different theory will be required.

Application to tunnelling through quantum dots

Although we derived our semiclassical expressions with low-temperature bridge-mediated electron transfer in mind, the fourth-order rate theory is quite general and should be valid for any set of PESs that fit the description above. We shall illustrate the new theory using a simple one-dimensional system which models three quantum dots in a row with nearest-neighbour coupling. Tunnelling plays a key role in describing various useful properties of quantum dots that are used in single-electron transistors [98] and as qubits for quantum computers [99].

The Born–Oppenheimer formalism typically treats nuclear dynamics moving on a PES or PESs defined by integrating out the electronic degrees of freedom. However, for a simple description of quantum dots, one typically treats a single electron as the dynamical variable and constructs potential-energy surfaces which describe its interaction with the bulk semiconductor. The mathematical formalism is otherwise identical. For small diabatic couplings, the adiabatic PESs will approach each other near the crossing points, thus invalidating the Born–Oppenheimer approximation and requiring a non-adiabatic rate theory.

Bound electron states of a quantum dot are commonly modelled as harmonic potentials [100] and so we define where denotes the coordinate of the single electron. We choose parameters that correspond to a self-assembled InAs quantum dot [101] with effective electron mass and vibrational energy spacing , which for simplicity is chosen to be the same for all states. The bridge and product minima are located at and , respectively, from the reactant, which corresponds to a typical length scale for quantum dots. The bias for the bridge state is set to while the bias of the product state, , is varied from 0 to 0.88 eV. Finally, we choose a temperature of 125 K for all our calculations and assume that the couplings and are independent of position.

For such a system, the bridge-mediated fourth-order rate describes a process called elastic cotunnelling [100,102], in which an electron tunnels from a source () to a drain ( via the non-resonant virtual states of an intermediate quantum dot (). These virtual states are higher in energy than the initial and final state, which themselves must be equal in energy for the reaction to take place. We can see that this exact condition is what is imposed by the single Dirac delta function in equation (4.4). At higher temperatures, a sequential process may become dominant, for which the instanton theory derived in this work is not applicable.

For the quantum-mechanical calculations, we use equation (4.4) and smear the delta function by replacing it with . This accounts for the physical mechanism of dissipation in the wells and allows an incoherent rate constant to be defined in a one-dimensional system which would otherwise display coherent dynamics. The substitution is equivalent to truncating the integral of the correlation function to the range . For simplicity, quantum-mechanical results were computed only for biases at which the reactant and product vibrational states matched exactly, for which the sinc function takes a value of , whereas it is 0 for any channel which does not exactly conserve energy.

The quantum-mechanical calculations required a relatively large number of bridge states to converge because the overlap terms oscillate as functions of the bridge level, as shown in figure 5. The small magnitude of the values also makes computing the sum numerically challenging. For this particular model system, standard double-precision floating-point numbers were found to be adequate, but the nature of the sum means that one may require higher precision arithmetic in some cases.

Figure 5

(a) The oscillations of the overlap between the reactant state and product state through the bridge state , defined as , as a function of the bridge state . Lines connecting integer values of are used only to ease visualization. (b) The exact quantum-mechanical rate as a function of the maximum bridge level, , included in the sum. The SCI rate is shown as a horizontal line as the instanton formalism does not require a sum over bridge states. All calculations were done for the model system with a bias of . (Online version in colour.)

For a particle in a harmonic potential, analytic expressions for the imaginary-time actions are known [6]. We use these results to evaluate the fourth-order instanton rates for our model system according to equation (4.10). According to the standard procedure, we use the semiclassical reactant partition function , which happens to be exact for harmonic systems.

The instanton for the system with a moderate bias of is shown in figure 4a. Its energy lies well below the energy of any bridge state, suggesting that the largest contributions to the tunnelling rate are from the lower vibrational states of . This is in agreement with the quantum-mechanical calculation, which indicates that the largest contribution comes from the vibrational ground state , followed by the first excited state .

The calculated rates are plotted as a function of the bias, , in figure 6a. Because the rate constant is proportional to the current, these results are related to the current–voltage plots often studied in molecular electronics [103]. We see that the SCI and exact rates generally agree very well with each other, are within about 20% at low biases and differ by a factor of approximately 2 at the largest bias considered. This point corresponds to an instanton with an infinitesimally small trajectory on and thus represents the limit of the method’s applicability. At stronger biases than this, a different SCI method would need to be derived, as the energies at the RB and BP hopping points ( and ) would no longer be below . The main cause of the error can be traced back to the steepest-descent integration over and in equation (4.6). As illustrated in figure 6b, this overestimates the integral, particularly for high biases when the instanton has small values of .

Figure 6

(a) Coupling-independent rates as a function of the bias, , using our new fourth-order instanton theory (equation (4.10)), the exact quantum-mechanical expression (equation (4.4)), the classical approximation (equation (4.14)), and the superexchange-mediated effective coupling approximation (equation (2.4a)) with effective coupling (equation (4.12)). (b) A plot showing how a steepest-descent integration over overestimates the rate as a result of erroneously including the region . The two-time correlation function (defined via a steepest-descent integral of equation (4.7) over ) with varied and fixed at the value at the stationary point is shown with a solid line and the corresponding Gaussian used by the steepest-descent approximation with a dashed line of the same colour. Each pair of lines of the same colour is scaled by its value at the stationary point. (Online version in colour.)

Comparison with superexchange-mediated effective-coupling approximation

A commonly used approximation to the fourth-order rate can be made by assuming that, relative to the large gap between reactant and bridge, one can ignore the vibrational energy-level structure of the bridge. All bridge states are thus assigned the same energy and the sums over the bridge levels in equation (4.4) can then be eliminated using the resolution of the identity. This gives an approach known as the superexchange-mediated effective coupling (SEC) approximation [49] and has a form equivalent to the second-order rate from to (equation (2.4a)) but with an effective coupling , that approximately accounts for the influence of the bridge where and are the zero-point energies of the reactant and bridge. Just as for the exact fourth-order rate, we evaluate the expressions using a smeared delta function and a quantum-mechanical partition function as in §4(a).

From figure 6a, we observe that the SEC approximation significantly underestimates the rates when compared to the exact results and is much less accurate than the fourth-order instanton. This behaviour can be easily explained within the instanton formalism, which could in principle be applied on top of the SEC approximation. In this case, the dominant tunnelling pathway would travel only on and states, as in the golden-rule case. The total abbreviated action (equation (2.16c)) along the SEC pathway would be , which is necessarily greater than the fourth-order abbreviated action for a given energy, , as the paths travelling on the bridge have a lower potential energy than they would in the reactant or product states. By the same argument, it can also be shown that . It follows that the action must be greater than and that therefore the SEC approximation predicts that tunnelling occurs at a lower energy than necessary, with a larger action than in the fourth-order instanton theory, typically leading to an underestimation of the rate.

Note that there are various alternative options for how to define the effective coupling. In particular, one could enforce detailed balance using the denominator or instead [49]. However, it is clear that any different choice simply shifts the whole curve, and although this may improve the results overall, it cannot be used to give uniformly accurate rates for both large and small biases.

Comparison with classical sequential mechanism

The expression for the classical rate of a one-dimensional bridge-mediated transfer in the limit of small couplings is [12] where is the classical partition function for the reactant and is the higher of the reactant/bridge or product/bridge intersections, i.e. . The Landau–Zener [44,45] hopping probability in the non-adiabatic limit is (which includes a factor of 2 to account for the second crossing of nonreactive paths [12]), where is the difference in the slopes and is the classical velocity at the reactant/bridge crossing point, . is defined analogously for the product/bridge intersection.

Assuming , the integral can be evaluated analytically to give where is the zeroth modified Bessel function of the second kind. Figure 6a shows how the classical rate varies with respect to the bias for the model system and we see that it fails to predict the trend of the rate with respect to the product bias. Note also that the expression in equation (4.14) spuriously tends to infinity when the crossing point energies and are equal.

Using the asymptotic relation as , we can gain some useful insights into the mechanism underlying the classical rate. When the quantity is large (i.e. for low temperatures and strong exothermic bias), the rate tends to an expression with the exponential term , and we identify as the activation energy. In fact, the expression reduces in this limit to the classical second-order rate (equation (2.23)) (for a reaction from to ) with a classical reactant partition function and an effective coupling given by At the relatively low temperature of our problem, for which , this limit is valid across all but the smallest exothermic biases considered. As for the SEC approximation, the classical sequential mechanism leads to a rate expression which can be written as the second-order rate with an effective coupling. However, in contrast to the SEC approximation, the effective second-order rate corresponds to a transfer from to instead of to and thus has an activation energy which is independent of exothermic bias. This explains why the predictions are contrary to the quantum-mechanical benchmark.

Conclusion

In this paper, we have explained how SCI theory can be generalized to study tunnelling effects in certain non-adiabatic chemical reactions. In particular, we have reviewed its application in the golden-rule limit, derived new closed-form expressions in separable cases and demonstrated that it can give accurate results where simpler methods based on global linear approximations fail. It will be interesting to see if the predictions of instanton theory will be significantly different from those presented in the literature and whether better agreement with experiment can be achieved in this way.

We have presented the working equations to evaluate golden-rule instanton theory within the ring-polymer formalism. In particular, we demonstrate how to treat translational and rotational degrees of freedom in a manner applicable to molecular systems in the gas phase. These methods are currently being applied to various molecular systems in our group and results will be published in forthcoming work.

Finally, we have generalized instanton theory beyond the golden-rule limit to treat the rate of a fourth-order reaction proceeding from reactants to products via a bridge. Unlike expressions based on the classical sequential mechanism, instanton theory was able to correctly predict the dependence of the rate on the bias using a deep-tunnelling pathway under the bridge. It gave good quantitative agreement with benchmark quantum-mechanical calculations and was significantly more accurate than the superexchange-mediated effective coupling approximation. It would be trivial to extend the approach for cases with multiple intermediate bridges. However, further extensions will be necessary to treat systems in the inverted regime or with high-energy bridges.

The nature of the systems for which this fourth-order instanton theory is applicable requires long tunnelling pathways. These will only exist at low temperature and for light tunnelling particles. For this reason, we chose to demonstrate the method using tunnelling of an electron through a row of quantum dots. It is clear that an accurate treatment of deep tunnelling is essential for this problem, as the classical mechanism fails to predict the correct trends.

In this work, we have applied the new fourth-order instanton approach only to a one-dimensional harmonic model system for which the propagators are known analytically. The ring-polymer approach used successfully for the second-order golden-rule case can, however, be directly generalized to the fourth-order rate, which would make this method applicable to complex multidimensional anharmonic systems.