Fermi’s golden rule defines the transition rate between weakly coupled states and can thus be used to describe a multitude of molecular processes including electron-transfer reactions and light-matter interaction. However, it can only be calculated if the wave functions of all internal states are known, which is typically not the case in molecular systems. Marcus theory provides a closed-form expression for the rate constant, which is a classical limit of the golden rule, and indicates the existence of a normal regime and an inverted regime. Semiclassical instanton theory presents a more accurate approximation to the golden-rule rate including nuclear quantum effects such as tunneling, which has so far been applicable to complex anharmonic systems in the normal regime only. In this paper, we extend the instanton method to the inverted regime and study the properties of the periodic orbit, which describes the tunneling mechanism via two imaginary-time trajectories, one of which now travels in negative imaginary time. It is known that tunneling is particularly prevalent in the inverted regime, even at room temperature, and thus, this method is expected to be useful in studying a wide range of molecular transitions occurring in this regime.

Describing chemical reactions that take place on more than one electronic potential-energy surface poses one of the primary open challenges in the field of chemical reaction dynamics.1,2 These processes are relevant to many phenomena, which we encounter not only in different disciplines of science but also in our everyday life. Ranging from redox reactions to photosynthesis, harvesting light in solar cells, molecular switches, and many more, the most fundamental step of these processes is a nonadiabatic transition from one electronic state to another, leading to a breakdown of the Born–Oppenheimer approximation.3 

Due to the great interest in these phenomena, the study of nonadiabatic transitions is an important topic for research. Hence, there exists a plethora of different algorithms to simulate nonadiabatic dynamics,3–8 from computationally expensive but accurate methods based on wave-function propagation9–12 to quantum-classical mapping approaches13–18 and heuristically motivated, pragmatic methods such as trajectory surface hopping.19–21 

Simulating the direct dynamics of a chemical reaction, however, is not usually a practical way to obtain information about the reaction rate because the typical time scales of chemical reactions are long. Instead, a nonadiabatic extension of transition-state theory (TST) is required.22 

The thermal rate constant for the transition from the reactant electronic state, with internal energy levels E0λ and a partition function Z0=λeβE0λ, to the product electronic state, with internal energy levels E1ν, can be found by applying perturbation theory. The result to lowest order in the coupling Δλν between these states is given by the famous golden-rule formula23,24 generalized for an initial thermal distribution,25,26
kQM=2πλeβE0λZ0ν|Δλν|2δ(E0λE1ν),
(1)
whose name (given by Fermi) indicates its overwhelming relevance and applicability in a multitude of different fields, including nuclear physics, light-matter interactions, and nonadiabatic transitions.27–31 One important example of the latter is the nonadiabatic transfer of an electron from a donor to an acceptor.32 

Marcus was awarded with the Nobel Prize in chemistry in 199232 for his work on electron-transfer rate theory.33,34 One of the great triumphs of his theory was the prediction of an inverted regime, in which the rates decrease despite increasing thermodynamic driving force, which was later confirmed by experiment.35 Because of its simplicity and practicability, Marcus theory remains the most commonly applied approach to describe electron-transfer reactions.36–42 

However, there are a number of approximations inherent in Marcus theory,43 which includes the assumption of parabolic free-energy curves along the reaction coordinate. It also employs classical statistics and cannot, therefore, capture nuclear quantum effects such as zero-point energy and quantum tunneling, the neglect of which could lead to deviations from the exact rate of several orders of magnitude, especially at low temperatures. The inclusion of these effects in novel nonadiabatic rate theories that can be applied to molecular systems without making the parabolic approximation is, therefore, a major objective.1 

In particular, it has been predicted that quantum tunneling effects can substantially speed up the rate in the inverted regime.44 This was confirmed experimentally in Ref. 45, in which reaction rates were found to be up to eight orders of magnitude larger than were predicted by Marcus theory but which could be explained by including a quantum-mechanical treatment of the vibrational modes.46 Early approaches such as these for including quantum statistics into Marcus theory were, however, only possible by restricting the system to simplistic models such as the spin-boson model.47 On the other hand, classical golden-rule transition-state theory43,48,49 has no restriction on the complexity of the system but does not take account of quantum nuclear effects.

A domain of methods that proved to be particularly successful in describing nuclear quantum effects in more general systems is based on Feynman’s path-integral formulation of quantum mechanics.50 For instance, Fermi’s golden rule has been recast into a semiclassical instanton (SCI) formulation,51 which does not require detailed knowledge about the internal states and can, therefore, be applied to complex systems. However, there is a problem with these methods in the inverted regime because the imaginary-time propagator, on which most of these methods rely, diverges. Hence, many previous semiclassical51–58 and imaginary-time path-integral approaches59–61 did not tackle the inverted regime.

One approach for extending these methods to treat the inverted regime was suggested by Lawrence and Manolopoulos,62 who analytically continued Wolynes’ rate expression,59 which was derived from a cumulant expansion of the correlation function, into the inverted regime. The rate is then obtained by extrapolating the path-integral data collected in the normal regime into the inverted regime. While their methodology appears to work very well at predicting rates and could, in principle, be directly applied to extend the instanton rate in a similar way, the mechanistic view may be lost by this approach, as the rate is not extracted directly from a simulation of the system in the inverted regime.

The electron-transfer rate in the inverted regime cannot be tackled directly by standard ring-polymer molecular dynamics63 where the transferred electron is treated as an explicit particle,64 as can be explained by a semiclassical analysis.65 However, more recent modifications of ring-polymer molecular dynamics66–68 and golden-rule quantum transition-state theory69,70 have started to address this problem, but these still lack the rigor, simplicity of implementation, and mechanistic insight which one can obtain from instanton theory.

In this paper, we propose an extension of the semiclassical instanton method51–53 for the inverted regime in the golden-rule limit. Building on the successful approach of Ref. 62, we analytically continue the rate expression derived for the normal regime to obtain a method applicable to the inverted regime. However, in contrast to Ref. 62, the instanton rate formula we obtain by this approach is a one-shot method, which requires no extrapolation of results collected in a different regime. We show excellent agreement with the exact golden-rule rates for the model systems studied. At the same time, it gives direct mechanistic insights, as it automatically locates the dominant tunneling pathway for the reaction under investigation, which is equivalent to predicting the mechanism. It can, therefore, be used to shed light on the role of quantum nuclear tunneling in electron-transfer reactions, which is expected to be of substantial importance, especially in the inverted regime.

The instanton approach can be implemented using a ring-polymer discretization, in which the only change necessary to make the algorithm applicable for the inverted regime is a slight variation in the optimization scheme, which turns out to be just as reliable and effective as the optimization in the normal regime. Hence, the method is conceptually ideally suited for use in conjunction with ab initio potentials and thereby for realistic simulations of molecular systems, as has been demonstrated for the standard ring-polymer instanton approach.71–77 

In this section, we summarize our previous derivation of semiclassical instanton theory for the golden-rule rate in the normal regime.51 This follows a similar approach to our derivation of instanton theory on a single Born–Oppenheimer potential.78,79

We consider a general multidimensional system with two electronic states, each with a nuclear Hamiltonian of the form
Ĥn=j=1Dp^j22m+Vn(x^),
(2)
where n ∈ {0, 1} is the electronic-state index and x = (x1, …, xD) are the Cartesian coordinates of D nuclear degrees of freedom. These nuclei move on the potential-energy surface Vn(x) with conjugate momenta p^j. Without loss of generality, the nuclear degrees of freedom have been mass-weighted such that each has the same mass, m. The electronic states |n⟩ are coupled by Δ(x^) to give the total Hamiltonian in the diabatic representation,43,
Ĥ=Ĥ0|00|+Ĥ1|11|+Δ(x^)|01|+|10|.
(3)
We shall take the diabatic coupling to be constant, Δ(x^)=Δ, but note that the extension of the theory to coordinate-dependent couplings is fairly straightforward. Furthermore, we assume the coupling to be very weak, i.e., Δ → 0, known as the golden-rule limit, which is typically the case in electron-transfer reactions.43 

The quantum-mechanical (QM) rate is then given by the golden-rule expression [Eq. (1)], which is valid both in the normal and inverted regimes. An illustration of a typical electron-transfer system is given in Fig. 1, showing the internal states involved at a particular energy E=E0λ=E1ν for which nuclear tunneling plays a role. However, in order to calculate the rate in this way, the internal states of both the reactant and product would be required, which are typically not known and cannot be computed for a complex system.

FIG. 1.

Illustration of a typical electron-transfer system with diabatic potentials, Vn(x), in color and the corresponding adiabatic potentials (calculated with a small value of Δ) depicted by the black dashed lines in the normal regime (a) and inverted regime (b). The nuclear wave functions ψ0λ(x) and ψ1ν(x) corresponding to eigenstates of energy E=E0λ=E1ν of the two diabatic-state Hamiltonians are shown in red and blue with an arbitrary normalization. The overlap of the wave functions is shown by the purple dotted line both in the main figures and the insets which enlarge the area framed by the gray boxes. The coupling Δλν is proportional to the integral over this function. The diabatic crossing point, x, and the corresponding potential energy, V, are marked on the axes.

FIG. 1.

Illustration of a typical electron-transfer system with diabatic potentials, Vn(x), in color and the corresponding adiabatic potentials (calculated with a small value of Δ) depicted by the black dashed lines in the normal regime (a) and inverted regime (b). The nuclear wave functions ψ0λ(x) and ψ1ν(x) corresponding to eigenstates of energy E=E0λ=E1ν of the two diabatic-state Hamiltonians are shown in red and blue with an arbitrary normalization. The overlap of the wave functions is shown by the purple dotted line both in the main figures and the insets which enlarge the area framed by the gray boxes. The coupling Δλν is proportional to the integral over this function. The diabatic crossing point, x, and the corresponding potential energy, V, are marked on the axes.

Close modal

It is clear, however, from Fig. 1(a) that, at least in the normal regime, the rate expression will be dominated by the behavior of the wave functions near the crossing point and that, therefore, it should be possible to find a good approximation relying only on local knowledge of the system close to this point. This is not immediately clear for the inverted regime shown in Fig. 1(b), although we will nonetheless be able to obtain a simple instanton expression for the rate in Sec. III.

The purpose of the semiclassical instanton approach is to obtain a good approximation to the golden-rule rate without detailed knowledge of the internal states. Therefore, instead of using the expression in Eq. (1), we employ the alternative exact definition of the quantum rate,43,48,80
kQMZ0=Δ22c(τ+it)dt,
(4)
where the flux correlation function is
c(τ+it)=TreĤ0(βτit)/eĤ1(τ+it)/
(5)
and the reactant partition function is Z0=TreβĤ0. Note that in order to write the expression in this form, it is necessary to assume that the energies of the internal states of both the reactant and product are bounded from below, i.e., there exists a finite-energy ground state of Ĥ0 and Ĥ1. We shall in addition assume that the energies are not bounded from above, which is the typical situation for molecular Hamiltonians.

In this section, we shall choose τ in the range 0 < τ < βℏ. Under these circumstances, it can be shown that c(iz) is an analytic function of z = t − iτ and as such the integral c(iz)dz is independent of the contour of integration, and hence, the rate is independent of the choice of τ, at least within its range of validity.

Expanding the trace in a coordinate-space representation gives
kQMZ0=Δ22K0(x,x,βτit)×K1(x,x,τ+it)dxdxdt,
(6)
where the imaginary-time quantum propagators, defined by
Kn(xi,xf,τn)=xf|eτnĤn/|xi,
(7)
describe the dynamics of the system evolving from the initial position xi to the final position xf in imaginary time τn according to the Hamiltonian Ĥn. Real-time dynamics can also be described by making the third argument complex.

Equation (6) is valid for systems in both the normal and inverted regimes and can be evaluated for model systems where the propagators are known analytically by numerical integration. However, because it is necessary to limit ourselves to the range 0 < τ < βℏ, as we will show, the semiclassical approximation described in Sec. II B can only be derived directly for the normal regime.

The instanton expression for the rate is obtained by first replacing the exact quantum propagators by semiclassical van-Vleck propagators81 generalized for imaginary-time arguments,79,82
Kn(xi,xf,τn)Cn(2π)DeSn/.
(8)
This expression is evaluated using the classical trajectory, xn(u), which travels from xn(0) = xi to xn(τn) = xf in imaginary time τn. This trajectory is found as the path that makes the Euclidean action, Sn, stationary, and the action is defined as
Sn=Sn(xi,xf,τn)=0τn12mẋn(u)2+Vn(xn(u))du,
(9)
where ẋn(u)=dxndu is the imaginary-time velocity. The prefactor of the semiclassical propagator is given by the determinant
Cn=2Snxixf.
(10)
Plugging this semiclassical propagator into Eq. (6) allows us to perform the integrals over the end-points x′, x″ and over time t, employing the method of steepest descent.83 This leads to the following expression for the golden-rule instanton rate:51,
kSCIZ0=2πΔ22C0C1ΣeS/,
(11)
where the total action is
S=S(x,x,τ)=S0(x,x,βτ)+S1(x,x,τ),
(12)
and the determinant arising from the steepest-descent integration is
Σ=2Sxx2Sxx2Sxτ2Sxx2Sxx2Sxτ2Sτx2Sτx2Sτ2.
(13)
All these expressions are evaluated at a set of values for x′, x″, and τ, which describes the stationary point of the action, defined by Sx=Sx=Sτ=0. If τ is chosen according to this prescription, it is not even necessary to deform the integration contour, as the saddle point appears at t = 0 with both x′ and x″ purely real. The minus sign before Σ in Eq. (11) arises naturally from the Cauchy–Riemann relations84 when re-expressing derivatives with respect to t as derivatives with respect to τ. If there is more than one stationary point of the action, the rate is given by a sum over each solution.70 For consistency, the reactant partition function, Z0, should also be evaluated within a semiclassical approximation.79 
If we instead take the steepest-descent integration over the end-points first and then separately over time, this leads to the alternative but equivalent expression85 
kSCIZ0=2πΔ22C0C1Cd2Sdτ212eS/,
(14)
where
C=2Sxx2Sxx2Sxx2Sxx.
(15)

Thus, the total action, S(x′, x″, τ) [Eq. (12)], is a sum of the actions of two imaginary-time classical trajectories, one for each electronic state. One trajectory travels on the reactant state from x′ to x″ in imaginary time τ0 = βℏτ and the other from x″ to x′ in imaginary time τ1 = τ on the product state. This forms a closed path of total imaginary time τ0 + τ1βℏ, known as the instanton.

Classical trajectories in imaginary time are described by ordinary classical mechanics but with an upside-down potential.82 They describe quantum tunneling by traveling through the classically forbidden region and typically “bounce” against the potential, which we define as an encounter with a turning point such that the momentum is instantaneously zero.51 

As has been shown in Ref. 51 for instantons in the normal regime, the fact that x′ and x″ are chosen as stationary points of S is to say that the imaginary-time momenta on each surface, pn(u)=mẋn(u), must have the same magnitude and point in the same direction at the end-points. Hence, the two classical trajectories join smoothly into each other. Furthermore, the restriction Sτ=0 ensures energy conservation, which implies that the instanton is a periodic orbit. Typically, the two end-points will be located in the same place, which we call the hopping point, x, and we can conclude that it must be located on the crossing seam of the two potentials, where V0(x) = V1(x). Further details about these statements will be given in Sec. III B.

All the steps in this derivation are asymptotic approximations, which become exact in the → 0 limit (with βℏ kept constant). This is, therefore, known as a semiclassical approximation. Semiclassical instanton theory gives the exact rate for a system of two linear potentials, and for more general systems in the limit of high temperature or heavy masses, it approaches a harmonic approximation to classical transition-state theory.51 In practice, the theory is applied using a ring-polymer discretization of the imaginary-time trajectories following the approach described in Ref. 52. This method has been previously used to study tunneling effects and the dependence of asymmetrical reorganization energies in a system-bath model in the normal regime.53 It has also been employed in quantifying the contributions from competing reaction mechanisms.70 

Here, we consider the classical, high-temperature limit (βℏ → 0) of a general curve crossing problem. The instanton for this system corresponds to two very short imaginary-time trajectories describing a periodic orbit, which is collapsed onto an infinitesimally short line.51 Note that unlike for instanton theory on a single Born–Oppenheimer potential,76,86 there is no cross-over temperature for golden-rule instanton theory.

For a one-dimensional system in this classical limit, the action of a single trajectory can be written in the simpler form
Sn(xi,xf,τn)=m2τnx2+Vn(x+)τn,
(16)
where x = xfxi and x+=12(xi+xf). The stationary points of the total action [Eq. (12)] can be found by searching first for the solution to Sx=0, which gives x = 0, and then for the solution to Sτ=0 evaluated at x = 0, which requires that the hopping point, x+ = x, obeys V0(x) = V1(x). Finally, Sx+=0 requires that
τ=βV0(x)V0(x)V1(x),
(17)
where ∇Vn(x) is the derivative of the potential Vn(x) with respect to x evaluated at x. These solutions give the simple interpretation that the transition-state for the reaction is located at the crossing point between the two diabatic potentials, x. Although the value of τ that makes S stationary does not have a clear interpretation within the classical theory, it plays an important role in defining the instanton. We shall, therefore, consider the behavior of τ in various regimes.

The common definition of the “inverted” regime is typically expressed in the context of Marcus theory by a system that has a larger driving force than reorganization energy. A more general definition is that the different regimes are defined by the slope of the potentials at the crossing point; in the inverted regime, the gradients have the same sign, whereas in the normal regime, the gradients have opposite signs. An alternative terminology for these two cases is “Landau–Zener” type or “nonadiabatic-tunneling” type.87 Note that the common definition is equivalent to the more general definition as long as the driving force is defined as ε=V0(xmin(0))V1(xmin(1)) and the (product) reorganization energy as Λ=V1(xmin(0))V1(xmin(1)), where xmin(n) is the minimum of Vn(x). In multidimensional systems, there is a crossing seam, and one would say that the scalar product of the two gradient vectors on this seam is positive only in the inverted regime. In fact, at the minimum-energy crossing point, which is the location of the hopping point in the classical limit, these gradient vectors are antiparallel in the normal regime51 but parallel in the inverted regime.

From Eq. (17), it can be seen that in the normal regime, where ∇V0(x) and ∇V1(x) have opposite signs, the value of τ that makes S stationary falls in the range 0 < τ < βℏ and is, therefore, always positive. In the activationless regime, where ∇V0(x) = 0, the situation changes to τ = 0. In this paper, we shall consider only one type of inverted regime, ∇V1(x)/∇V0(x) > 1, which is the typical case encountered in Marcus theory. Here, τ takes a negative value. Needless to say, all our results could be easily converted to describe a system in the alternative inverted regime, ∇V0(x)/∇V1(x) > 1. Equation (17) generalizes to multidimensional systems by projecting each gradient vector along the same arbitrary direction, and τ is thus seen to have the same behavior.

Therefore, in the context of the high-temperature limit of a general (anharmonic) system of two potential-energy surfaces that cross, we have shown that the sign of τ is different in the two regimes. This rule is also known to hold for situations involving quantum tunneling,48,62,88,89 which will be confirmed by the analysis later in this paper.

When substituting the action [Eq. (16)] evaluated at its stationary point into Eq. (11), one obtains the classical rate
kclTSTZ0cl=2πmβ2Δ2|V0(x)V1(x)|eβV,
(18)
where V = V0(x) = V1(x) and Z0cl is the classical limit of the reactant partition function. Note that this expression for the rate is equivalent to that derived from classical statistical mechanics employing the Landau–Zener hopping probability in the golden-rule limit.26,48 It can also be derived directly from classical golden-rule TST,43,48 which is proportional to eβV0δ(V0V1)dx, by noting that the integral can be performed easily for one-dimensional systems due to the constraint. For a spin-boson model, this classical expression reduces to Marcus theory.

Equation (18) is, in fact, valid, not just in the normal regime but also for the inverted regime, that is, whether ∇V0(x) and ∇V1(x) have the same sign or opposite signs, as long as they are not equal to each other. It is noteworthy that we can obtain valid classical formulas for the inverted regime from this approach even though in the derivation we assumed that 0 < τ < βℏ, which is only appropriate in the normal regime. In Sec. III, we shall also attempt to generalize the instanton method to the inverted regime in a similar way.

There are a number of differences between the normal and inverted regimes, some of which can be seen directly in Fig. 1. In the normal regime, the overlap between the wave functions, emphasized in the insets, has a shape similar to a Gaussian centered around the diabatic crossing point x. As instanton theory is derived by performing integrals by steepest descent, it requires the integrands to have an approximately Gaussian shape. It is, therefore, not surprising that instanton theory can successfully be applied to treat the normal regime.51 In the inverted regime, however, the shape of the overlap is not so localized and is oscillatory, raising the question if the semiclassical approximation is valid in this regime at all.

Another difference is that in the normal regime, the adiabatic ground-state potential has a cusp shape, which exhibits a distinct barrier separating reactants and products, as illustrated by the black dashed lines in Fig. 1(a). In fact, the instanton orbit on the ground-state adiabatic cusped potential can be used to find the golden-rule instanton in the normal regime.54 This procedure will not, however, work in the inverted regime for which no barrier on the ground-state adiabatic potential exists. In fact, the rate within the Born–Oppenheimer approximation is not even defined, as the ground-state adiabatic potential only has one minimum.

Probably the most significant difference is that in Sec. II C, we defined the inverted regime using the gradients on the crossing seam and found that the value of τ that makes S stationary becomes negative. With the definitions given above, this implies τ1 < 0 and τ0 > βℏ in such a way that τ0 + τ1βℏ still holds. This has important consequences for the implementation and interpretation of the instanton approach, which we discuss in this section.

The main complication for the derivation of instanton theory in the inverted regime is that the expression for the imaginary-time propagator K1 diverges for negative τ1. One cannot, therefore, write the rate in terms of the coordinate-space integral as in Eq. (6) using the appropriate value for τ, which would be necessary in order to carry out the steepest-descent integration. However, we will show that while this path-integral expression does diverge, Fermi’s golden rule remains well defined and can be approximated to good accuracy by an analytic continuation of the instanton rate formula.

Here, we study more carefully the cause of the divergence in the inverted regime. To do this, we shall investigate the correlation function written as a sum over the contributions of the eigenstates, ψ0λ, of the reactant Hamiltonian, Ĥ0, labeled by the quantum number λ and the eigenstates, ψ1ν, of the product Hamiltonian, Ĥ1, labeled by ν.90 The flux correlation function [Eq. (5)] expanded in the energy basis is thus given by
c̃(τ)=νgν(τ),
(19a)
gν(τ)=λ|θλν|2e(βτ)E0λ/τE1ν/,
(19b)
where θλν=ψ0λ(x)*ψ1ν(x)dx such that the coupling between states used in Eq. (1) is given by Δλν = Δθλν. The overlap integral is clearly bounded by 0 ≤ |θλν| ≤ 1, which follows from the Cauchy–Schwarz inequality assuming that the wave functions are square-integrable and normalized. Examples of these overlap integrals are depicted in Fig. 1, although note that in Eqs. (19), there is no restriction that the energies of the two internal states must be the same.

In order to discuss the convergence of the correlation function, let us first assume that we have a finite system for which the partition functions exist at any temperature, i.e., the sums νeτnEnν/ converge for any τn > 0. In both the normal and inverted regimes, we are interested in values of βℏτ > 0, and so it is clear from the comparison test that the sum over λ converges, making gν(τ) a well-defined quantity in either regime. By similar arguments, the sum over ν is also clearly convergent when τ is in the range 0 < τ < βℏ, which would be appropriate only for the normal regime.

Let us also consider the flux correlation function expanded in the position basis,
c(τ)=K0(x,x,βτ)K1(x,x,τ)dxdx,
(20)
where Kn(xi,xf,τn)=νψnν(xi)*ψnν(xf)eτnEnν/. We expect that for a physical system, |ψnν(x)| is bounded by a positive number, Ψmax, for all ν and x. Therefore, the absolute value of any term in the sum is necessarily less than Ψmax2eτnEnν/. The comparison test (also known as the “Weierstrass M test”)84 can again be invoked to prove that the propagators converge for τn > 0 for any values of x′ and x″, which is known as uniform convergence.84 Note, however, that K1(xi, xf, τ1) diverges in the inverted regime if we choose τ1 < 0, as is required to perform the steepest-descent integration.

Inserting the definition for θλν into c̃(τ) and the wave-function expansion for Kn into c(τ), it is clear that the two correlation functions are defined in a similar way, the only difference being that the sums and integrals are taken in a different order. Only for a uniformly convergent series can we interchange sums and integrals without affecting the result, and thus, it is possible to show that c̃(τ)=c(τ) only for 0 < τ < βℏ. Because c̃(iz) and c(iz) are analytic functions of z = t − iτ in the regions where they converge, this analysis can easily be extended to study the correlation functions at any value of t. This simply adds phases to each term, which does not change the convergence behavior and the fact that they are identical for 0 < τ < βℏ.

If we choose τ to be negative, however, as would be appropriate for the case of the inverted regime, the sum in K1 is no longer uniformly convergent and thus c(τ) diverges. Interestingly, we find that the correlation function c̃(τ) remains (at least in some cases) well defined. To demonstrate this, we take as an example the one-dimensional version of the spin-boson model, defined in Sec. V A, deep in the inverted regime with a driving force twice that of the reorganization energy. Using the value of τ that makes S stationary, we evaluate the terms in Eq. (19) in the eigenfunction bases of the two harmonic oscillators. In Fig. 2(a), we plot the contributions from each term in the series and demonstrate that the gν(τ) terms exhibit a distinct peak and falls off rapidly either side. This occurs because θλν exponentially decreases for states of widely different energies. Therefore, the correlation function c̃(τ) converges in this example and is analytic even for systems in the inverted regime where τ is negative.

FIG. 2.

(a) Terms contributing to the trace defined by Eq. (19), each of which corresponds to a product eigenstate with quantum number ν plotted against quantum number ν and corresponding energy E. (b) Plot of the real (blue solid line) and imaginary parts (orange dotted line) of the flux correlation function. In addition, the steepest-descent approximation to the flux correlation function is shown (black dashed line). All quantities are computed for the spin-boson model of Sec. V A with D = 1 and ε/Λ = 2 for τ ≈ −0.26βℏ, which is the stationary point of the action for this system.

FIG. 2.

(a) Terms contributing to the trace defined by Eq. (19), each of which corresponds to a product eigenstate with quantum number ν plotted against quantum number ν and corresponding energy E. (b) Plot of the real (blue solid line) and imaginary parts (orange dotted line) of the flux correlation function. In addition, the steepest-descent approximation to the flux correlation function is shown (black dashed line). All quantities are computed for the spin-boson model of Sec. V A with D = 1 and ε/Λ = 2 for τ ≈ −0.26βℏ, which is the stationary point of the action for this system.

Close modal

In Fig. 2(b), we show the dependence of the flux correlation function on real time, t, at a fixed negative value of τ, which demonstrates a simple form in almost perfect agreement with the Gaussian function employed by steepest descent. Hence, even though the overlap function used to calculate the golden-rule rate in the time-independent picture represented by Eq. (1) and illustrated in Fig. 1 is delocalized and oscillatory, the correlation function at the optimal τ can give a simple physical picture of the tunneling. The peak in Fig. 2(a) is centered around an energy, which will later turn out to be the instanton energy (defined in Sec. III B), and the correlation function in Fig. 2(b) is dominated by short-time behavior and still has an approximately Gaussian shape, justifying the use of steepest-descent integration.91 That this is possible is, of course, no surprise as the analytic continuation of the spin-boson formula and Wolynes theory also relies on this behavior.62 

In the normal regime, we know how to make good semiclassical approximations to c(τ) using instanton theory but have no simple approach based on c̃(τ), which is defined in terms of wave functions. Therefore, in order to formulate instanton theory in the inverted regime, we employ the mathematical process of analytic continuation.84 This allows us to evaluate an approximation to c(iz) for positive τ, which, because c(iz)=c̃(iz), must also be a good approximation to c̃(iz) in this regime. Because c̃(iz) is analytic across both regimes, this approximation will also be valid in the inverted regime. Accordingly, we propose to analytically continue the instanton method into the inverted regime and will employ the semiclassical instanton rate expression of Eq. (11), not just in the normal regime, where it was originally derived, but also for the inverted regime.

In effect, the method of Ref. 62 analytically continued the function c(τ) into the region with τ < 0 by fitting it numerically to a suitable functional form2 based on calculations in the normal regime and extrapolating to describe systems in the inverted regime. We build on this approach but go one step further to find a semiclassical instanton approach that is directly applicable in the inverted regime and requires no extrapolation. In the following, we analyze this new approach and show that it gives a valid approximation to Fermi’s golden-rule rate.

Through analytic continuation, we have a formula [Eq. (11)] for the golden-rule rate in the inverted regime based on an action S(x′, x″, τ). This should be evaluated at its stationary point, which defines the instanton and in this regime has a negative value for τ. In this section, we shall study the behavior of the instanton in the inverted regime and establish that it remains a periodic orbit, which travels through classically forbidden regions.

We start with the imaginary-time Euclidean action of a single trajectory, Sn(xi, xf, τn), defined by Eq. (9), and will be careful to ensure that all our formulas are valid for both positive and negative imaginary times, τn. The trajectory of interest, xn(u), starts at xn(0) = xi and travels to xn(τn) = xf in imaginary time τn. This trajectory has a conserved energy, En, because the Hamiltonian is time-independent. We can, therefore, add a zero under the integral
Sn(xi,xf,τn)=0τn12mẋn(u)2+Vn(xn(u))E(u)+Endu,
(21a)
where E(u)=12mẋn(u)2+Vn(xn(u)) is the instantaneous energy, which is constant (independent of u) and equal to En. Inserting this definition into Eq. (21a) leads to
Sn(xi,xf,τn)=0τnmẋn(u)2+Endu
(21b)
=xixfpndxn+Enτn,
(21c)
where dxn is an infinitesimal displacement vector pointing along the path in the direction from xi to xf. We call this the direction of particle flow. Our convention will be to define the imaginary-time momentum, pn=mdxdu, such that it points along the direction of change of position in a positive imaginary-time interval. Therefore, for a trajectory traveling from xi to xf in positive imaginary time τn > 0, the momentum, pn(u), will point along the path in the direction from xi to xf. However, for a trajectory traveling xi to xf in negative imaginary time τn < 0, the momentum, pn(u), will point in the opposite direction, i.e., along the path in the direction from xf to xi.

From these equations, we can determine that Sn(xi, xf, τn) ≡−Sn(xi, xf, −τn). This can be seen from Eq. (21b) or Eq. (9) as the integral changes sign when the integration range goes from 0 to a negative number. Alternatively, one can see that both terms in Eq. (21c) change sign when τn < 0 because for negative imaginary times, the momentum vector, pn, points in the opposite direction from the particle flow, i.e., it is antiparallel to dxn. In particular, if the zero potential is chosen below the instanton energy (for instance, at the reactant minimum), making En ≥ 0, then Sn will be positive when τn > 0 and negative when τn < 0. Therefore, whereas S0 remains positive just as in the normal regime, in the inverted regime, the value of S1 becomes negative.

As the instanton corresponds to a stationary point of the total action, we need to know the derivatives of the individual actions of each trajectory with respect to the initial and final end-points as well as with respect to imaginary time. These can be found by taking derivatives of Eq. (21c) to give
Sn(xi,xf,τn)xi=pn(0),
(22a)
Sn(xi,xf,τn)xf=pn(τn),
(22b)
Sn(xi,xf,τn)τn=En.
(22c)
Note that the derivative with respect to the initial point has a different sign from that with respect to the final point.81 
By utilizing these relations in the definition of the total action of the closed orbit [Eq. (12)], we arrive at the conditions
Sx=p0+p1,
(23a)
Sx=p0p1,
(23b)
Sτ=E0+E1,
(23c)
where pn is the momentum of the trajectory xn at the end marked x′ (and likewise for double primes), i.e., p0=mẋ0(0), p0=mẋ0(τ0), p1=mẋ1(τ1), and p1=mẋ1(0). All of the derivatives in Eqs. (23) must simultaneously vanish at the instanton configuration, which effectively imposes energy and momentum conservation at the intersection of the trajectories, p0 = p1, p0 = p1, and E0 = E1. The simplest solution to these equations (and typically the only one) is found when both x′ and x″ are located at the same coordinate, which we call the hopping point, x. The first two conditions require that at the hopping point, the momenta p0 and p1 must be vectors with the same direction and the same magnitude. Because the third condition requires the energies of the trajectories to match, E0 = E1, and to be conserved along the path, the potentials at the hopping point must be identical as well. We can thus conclude that the hopping point is located somewhere on the crossing seam between the two potential-energy surfaces such that V0(x) = V1(x). These findings are equivalent to those found in previous work limited to the normal regime,51 but we have now shown that they also hold in the inverted regime. However, there is nonetheless a fundamental difference in the inverted regime when we study the paths followed by the trajectories due to the fact that one path travels in negative imaginary time.

The imaginary-time classical trajectories, xn(u), that start and end at the same point x but travel in a nonzero amount of time τn, whether positive or negative, will typically bounce against the potential. This happens halfway along at the turning point, xnb=xn(τn/2), at which the kinetic energy is zero and the total energy En=Vn(xnb). These considerations, along with the conditions in Eq. (23), give us the picture shown in Fig. 3.

FIG. 3.

Visualization of the two imaginary-time trajectories forming instantons in (a) the normal and (b) the inverted regime at energy E = E0 = E1. The reactant trajectory, x0, is given in blue, and the product trajectory, x1, is given in red. The black arrows indicate the direction of particle flow from the initial point to the final point of each trajectory. The steepest-descent integration of positions is taken about the crossing point x′ = x″ = x at which V0(x) = V1(x) = V. In the normal regime, the trajectories bounce on either side of the crossing point, i.e., x0b<x<x1b, whereas in the inverted regime, both trajectories bounce on the same side of the crossing seam, i.e., x<x1b<x0b.

FIG. 3.

Visualization of the two imaginary-time trajectories forming instantons in (a) the normal and (b) the inverted regime at energy E = E0 = E1. The reactant trajectory, x0, is given in blue, and the product trajectory, x1, is given in red. The black arrows indicate the direction of particle flow from the initial point to the final point of each trajectory. The steepest-descent integration of positions is taken about the crossing point x′ = x″ = x at which V0(x) = V1(x) = V. In the normal regime, the trajectories bounce on either side of the crossing point, i.e., x0b<x<x1b, whereas in the inverted regime, both trajectories bounce on the same side of the crossing seam, i.e., x<x1b<x0b.

Close modal

In this way, we have discovered the form of the instanton appropriate for describing Fermi’s golden-rule in the inverted regime. We did this purely through a consideration of the stationary point of S, defined by Eq. (12). Like the instanton in the normal regime, it is localized in the vicinity of the crossing point and bears little resemblance to the delocalized wave-function overlap shown in Fig. 1(b). However, this instanton has a number of important differences when compared with that in the normal regime, as can be seen from the plots in Fig. 4 obtained by joining the two trajectories together to make the instanton orbit.

FIG. 4.

The two trajectories x0 (blue) and x1 (red) that form the instanton are plotted as a function of imaginary time, u, in (a) the normal and (b) the inverted regime. In both cases, the periodicity of the full instanton is βℏ and three cycles are shown. In the normal regime, the trajectories travel in positive imaginary time and bounce on opposite sides of x. In contrast, in the inverted regime, both trajectories bounce on the same side and in order to have a continuous path, x1 must travel backwards in imaginary time. The arrows indicate the direction of particle flow followed by trajectories from their initial point to their final point.

FIG. 4.

The two trajectories x0 (blue) and x1 (red) that form the instanton are plotted as a function of imaginary time, u, in (a) the normal and (b) the inverted regime. In both cases, the periodicity of the full instanton is βℏ and three cycles are shown. In the normal regime, the trajectories travel in positive imaginary time and bounce on opposite sides of x. In contrast, in the inverted regime, both trajectories bounce on the same side and in order to have a continuous path, x1 must travel backwards in imaginary time. The arrows indicate the direction of particle flow followed by trajectories from their initial point to their final point.

Close modal

In the normal regime, the dynamics are periodic with a periodicity of βℏ. The motion can be described as a particle that travels on the reactant state for an imaginary time τ0 and then suddenly turns into a particle on the product state with the same energy and momentum where it travels for an imaginary time τ1 before turning back into a reactant-state particle.

In the inverted regime, the instanton formed by joining the two trajectories is not single valued at certain times. However, we will still talk about it as an orbit because it remains periodic in imaginary time with periodicity βℏ and has continuous position and momentum vectors as well as conserving energy along its whole path. In particular, the energy and momentum are conserved at the two hopping points even though the path itself takes a sharp change of direction when the particle starts to travel in negative imaginary time. There is a similarity with these pictures and those used to explain the scattering of particles and antiparticles according to the Feynman–Stückelberg interpretation.92 The dynamics of antiparticles are equivalent to those of ordinary particles except that time is reversed so that they are found at their final point at an earlier time than their initial point.93 Reading Fig. 4(b) from left to right, one sees a single particle traveling on the reactant electronic state. At imaginary time u = 0, a new particle/antiparticle pair is created at the hopping point, while the old particle coexists at a different location. The new particle also travels on the reactant state but the antiparticle on the product state. At u = |τ|, the antiparticle annihilates with the original reactant particle at the hopping point, x, leaving only the new reactant particle, which continues in the role of the original particle when the cycle repeats.

The stationary point that corresponds to the instanton can be classified by an index, which is defined as the number of negative eigenvalues of the second-derivative matrix of the action. Knowing this index will be helpful not only for designing an algorithm for finding the instanton but also in order to compute the rate in the inverted regime, for which we need to compute determinants of second derivatives of the action. In the normal regime, C0 and C1 are always positive because both trajectories are minimum-action paths and Σ is negative because the instanton is a single-index saddle point in the space (x′, x″, τ). On the other hand, in the inverted regime, the trajectory x1 is a maximum-action path, and thus, all eigenvalues of the matrix 2S1xx are negative. C1, which is the determinant of this matrix, may thus be positive or negative depending on whether there are an even or an odd number of nuclear degrees of freedom, D. To find the signs of the second derivatives of the total action, we turn to the classical limit, which was studied in Ref. 51. From Eq. (67) in that paper, one can see that Σ has D + 1 negative eigenvalues and D positive eigenvalues in the inverted regime. As the instanton moves smoothly from the high-temperature classical limit to a low-temperature tunneling mechanism, there is no reason why the number of negative eigenvalues of Σ should change, so this result holds also in the general case. Therefore, Σ will always have the opposite sign from C1, ensuring that the square root in the prefactor of Eq. (11) remains real. Hence, the same instanton rate expression can be uniformly applied across both the normal and inverted regimes. Finally, we conclude that the instanton is a (D + 1)-index saddle point of S(x′, x″, τ), although due to the fact that x1 is a maximum-action path, it will have an even higher index in the space of paths as explained in Sec. IV.

Up to this point, we have employed the Lagrangian formulation of classical mechanics to define the imaginary-time trajectories. An alternative approach is provided by the Hamilton–Jacobi formulation in which the action is a function of energy rather than time.81 This leads to further insight into the behavior of the instanton in the inverted regime.

To derive the energy-dependent action, we start from Eq. (21c) and write S1 in a slightly different way to give
S1(xi,xf,τ1)=xixfp¯1dx+E1τ1.
(24)
By introducing the antiparticle momentum p¯1=p1, we align the direction of the momentum with the particle flow and thereby make the integral strictly positive. The magnitude of the imaginary-time momentum is pn(x,En)=2m[Vn(x)En], which is always real and positive in the classically forbidden region where the instanton exists. These definitions enable us to make the transition from the Lagrangian to the Hamilton–Jacobi formalism by defining the abbreviated actions as
Wn(xi,xf,En)=pn(x(s),En)ds,
(25)
where we have introduced the line element ds, which defines a metric in the configuration space such that (ds)2 = dx · dx. This gives a line integral along the path, defined in such a way that ∫ds is the path length. Note that Wn is purely positive and is independent of the sign of τn or the direction of the trajectories, which has already been accounted for by the sign in Eq. (24). This definition is, therefore, symmetric to an exchange of end-points, i.e., Wn(xi, xf, En) = Wn(xf, xi, En). The relations Sn = ±Wn + Enτn, where the sign is chosen to match the sign of τn, can thus be viewed as Legendre transforms obeying the conditions Snτn=En and WnEn=|τn|. It also follows that Snxi=±Wnxi and Snxf=±Wnxf.

It is well known that classical trajectories can be defined either as paths that make Sn stationary (known as Hamilton’s principle) or paths that make Wn stationary (known as Maupertuis’ principle). Typically, classical trajectories are minimum-action paths,94 and in the normal regime, both S0 and S1 are minima with respect to variation of the path. However, in the inverted regime, the product trajectory will be a maximum of S1. Trajectories that bounce once have a conjugate point81 and so the associated Wn are single-index saddle points with respect to variation of the path.51 There is nothing different in the definition of W1 in the inverted regime, so the product trajectory is also a single-index saddle point.

If we define WW(x′, x″, E) appropriately in the normal and inverted regimes, the total action [Eq. (12)] of the instanton can also be written as
S(x,x,τ)=W(x,x,E)+βE,
(26)
where either E = E0 = E1 is a function of (x′, x″, τ) via En=Snτn or β is a function of (x′, x″, E) via β=WE.
In order to clarify the definitions in each regime, we give an overview over the most important equations
NormalregimeInvertedregime0<τ0<β;0<τ1<βτ0>β;τ1<0S=S0(x,x,τ0)+S1(x,x,τ1)S=S0(x,x,τ0)+S1(x,x,τ1)S0(xi,xf,τ0)=W0(xi,xf,E0)+E0τ0S0(xi,xf,τ0)=W0(xi,xf,E0)+E0τ0S1(xi,xf,τ1)=W1(xi,xf,E1)+E1τ1S1(xi,xf,τ1)=W1(xi,xf,E1)+E1τ1W=W0(x,x,E)+W1(x,x,E)W=W0(x,x,E)W1(x,x,E),
where the definition of W follows in each case from the required relationship in Eq. (26). Differentiating Eq. (26) using the chain rule, we find Sx=Wx and Sx=Wx, which shows that the instanton, which is a stationary point of S, could also be defined as a stationary point of W with respect to x′ and x″. Either the corresponding temperature can be found for a given E using β=WE or E can be varied until this equation is solved for a given β.52 

Let us now check that these definitions are consistent with the one-dimensional schematic shown in Fig. 3. In the normal regime, it is clear that if we change either x′ or x″, we increase the path length of one trajectory while simultaneously decreasing the path length of the other. This increases one of the abbreviated actions, Wn, and decreases the other. In the normal regime, the derivative of the total abbreviated action, W = W0 + W1, vanishes at the hopping point because here W0x=W0x=W1x=W1x=p, where p = p0(x, E) = p1(x, E). In the case of the instanton in the inverted regime, changing the positions of the terminal points x′ or x″ leads to either an elongation or contraction of both trajectories. In this case, the total abbreviated action has been defined as W = W0W1 and its derivative also vanishes at the hopping point because here W0x=W0x=W1x=W1x=p.

Furthermore, it is possible to show in this formulation that in the normal regime, 2Wxx=2Wxx=2m[V0(x)V1(x)]/p>0, whereas in the inverted regime, 2Wxx=2Wxx=2m[V0(x)V1(x)]/p<0, and in both cases, 2Wxx=0. Therefore, the instanton is defined by a minimum of W(x′, x″) in the normal regime, but a maximum in the inverted regime.

Finally, we wish to check that the analytically continued instanton formula gives a reasonable physical model of tunneling in the inverted regime. Neglecting the less-important prefactors, the rate is proportional to eS/ = eβE eW/, where we have used Eq. (26) to separate the exponential into two terms. In a similar way to the standard instanton approach,76,86 we can interpret the first term as the probability of the system acquiring energy E from thermal fluctuations and the second term as being the probability of tunneling at this energy.

In order to make this interpretation, we require that W is positive, ensuring that the tunneling probability is bounded by 1. Noting that Wn ≥ 0, this is clearly the case in the normal regime, where W = W0 + W1. However, in the inverted regime, where W = W0W1, we will need to show that W0W1. This can be justified, at least for a system similar to the schematic in Fig. 3(b), for which x1 follows the part of the same path as x0, although it is expected to hold more generally. Using the fact that V0(x) ≥ V1(x) and thus p0(x, E0) ≥ p1(x, E1) for any point x along the instanton and noting also that the path length on the reactant state is longer, it is easy to see from the definition of Wn in Eq. (25) that W0W1 in this case.

As a corollary to finding that W ≥ 0, by taking the potential to be zero at the reactants, such that V0(xmin(0)), which ensures that E ≥ 0 for the instanton, it follows that S ≥ 0. This suggests that, at least within the approximation of neglecting the prefactors, the fastest rate will be found when S = 0, which is the activationless regime. There is a barrier to the reaction (relative to the reactant minimum) in both the normal and inverted regimes, which gives a positive action and thus a smaller rate. One thus recovers the Marcus turnover picture. However, because in the inverted regime, W is a difference rather than a sum of two contributions, it will typically be smaller than in the normal regime (for a given energy), leading to a larger tunneling effect. The optimal tunneling energy may also drop due to the longer imaginary-time taken by the reactant trajectory in the inverted regime. This increases S and, therefore, typically ensures that S is also larger in the inverted regime than the normal regime, again leading to a larger tunneling effect. This is certainly the case for the archetypal case of the system with linear crossed potentials. Here, the tunneling factor is known to be51  kSCI/kcl=exp[β32κ02κ12/24m(κ0κ1)2], which is seen to be larger in the case that the gradients κn have the same sign than if they have opposite signs for the same magnitudes.

For the standard instanton approach in the Born–Oppenheimer limit, it is known that the eW/ factor can be derived from one-dimensional WKB theory as the probability of reaction at a given energy, E.76,86 WKB theory has also been applied to the nonadiabatic transfer problem, and the result is known as the Zhu–Nakamura formula.87,95 In the diabatic limit, their result (written in our notation) for the inverted regime is e(W0W1)/, where the Wn values are calculated along paths from the crossing point to the turning points and back again. This is identical to the factor found here from analytic continuation of instanton theory and confirms that our new formula provides the correct physical description of the reaction according to the mechanistic interpretation suggested in Fig. 3.

Note that our formula can be applied to general multidimensional systems, whereas the Zhu–Nakamura formula is a one-dimensional theory. Therefore, only with instanton theory will it be possible to study the dominant tunneling pathways in multidimensional systems, which, in general, will involve all degrees of freedom and will typically not require that the x1 trajectory follows the part of the same path as x0.

In practical implementations of the instanton approach, we typically adopt a ring-polymer formulation.56,76,96–98 For golden-rule instanton calculations, we follow the approach suggested for the normal regime in Ref. 52 in which both paths are discretized into Nn equally spaced imaginary-time intervals of length ηn = τn/Nn. This same approach can be adapted for use in the inverted regime, where τ1 and, therefore, η1 will be negative, while N1, of course, remains positive. The resulting N = N0 + N1 beads describe the ring polymer as shown in Fig. 5.

FIG. 5.

Schematic showing the ring polymer corresponding to Eq. (27) for an example with N = 10 split into N0 = 6 and N1 = 4. There is no fundamental difference in the setup of the ring polymer between two regimes, but for clarity, we show the configurations into which the beads will automatically arrange themselves for (a) the normal and (b) the inverted regime. The beads are shown as circles colored blue if they feel the reactant potential, V0, and red if they feel the product potential, V1. Beads 0 and N0 feel an averaged potential. The springs between beads are represented by lines colored blue for an imaginary-time interval of η0 and red for an imaginary-time interval of η1. In the inverted regime, η1 is negative and, thus, the springs are repulsive.

FIG. 5.

Schematic showing the ring polymer corresponding to Eq. (27) for an example with N = 10 split into N0 = 6 and N1 = 4. There is no fundamental difference in the setup of the ring polymer between two regimes, but for clarity, we show the configurations into which the beads will automatically arrange themselves for (a) the normal and (b) the inverted regime. The beads are shown as circles colored blue if they feel the reactant potential, V0, and red if they feel the product potential, V1. Beads 0 and N0 feel an averaged potential. The springs between beads are represented by lines colored blue for an imaginary-time interval of η0 and red for an imaginary-time interval of η1. In the inverted regime, η1 is negative and, thus, the springs are repulsive.

Close modal
As previously discussed, the instanton corresponds to a stationary point of the action, which for a path described by a ring polymer is given by
SRP(x(0),,x(N1);τ)     =i=1N0mx(i)x(i1)22η0+i=1N01η0V0(x(i))        +i=N0+1Nmx(i)x(i1)22η1+i=N0+1N1η1V1(x(i))        +η0V0(x(0))+V0(x(N0))2+η1V1(x(N0))+V1(x(N))2,
(27)
where cyclic indices are implied such that x(0)x(N) in order to form a closed orbit. This is defined such that in the N0, N1 → ∞ limit, the value of the ring-polymer action at its stationary point is equal to the semiclassical instanton action, S.99 

As in the normal regime, we define the instanton as a stationary point of SRP with respect to the bead positions, {x(0), …, x(N−1)}, and τ simultaneously. Note that following our ring-polymer instanton approach for the normal regime,52 but in contrast to other path-integral approaches,59,62 here the ring-polymer spring constants or equivalently the imaginary-time intervals ηn do not have to be the same for the reactant and product beads. Hence, the optimization algorithm can vary τ independently of the N0:N1 split, which is chosen arbitrarily and kept fixed throughout the optimization. An intelligent choice can enhance convergence with respect to the number of beads but is not necessary to obtain a consistent and accurate rate prediction, as will be described in Sec. V. In order to further reduce the number of evaluations of the potential energy and its gradient, one could alternatively employ a half-ring polymer formalism in a similar way to the standard approach.52,96 Although this is valid for both the normal and inverted regimes, to simplify the explanations of our findings, we shall not take this extra step here.

In the normal regime, the instanton corresponds to a single-index saddle point. This is because the instanton is formed of two minimum-action paths, and S is a maximum only in the τ variable, i.e., d2Sdτ2<0. However, in the inverted regime, the stationary point of the action that we associate with our semiclassical instanton is not a single-index saddle point. In the inverted regime, we seek a saddle point with N0D positive and N1D + 1 negative eigenvalues. This can be understood by first considering the situation where the values of τ as well as x(0) and x(N0) are fixed at the hopping point. We would need to minimize the action with respect to the beads {x(1),,x(N01)} and maximize it with respect to the beads {x(N0+1),,x(N1)} in order to reproduce the instanton configuration. Note that this gives perfectly reasonable trajectories, as maximizing the action with respect to the second set of beads, which due to the fact that η1 < 0 have repulsive springs, is equivalent to minimization with attractive springs. Then, as explained in Sec. III B, the variation of the remaining points gives D negative and D positive eigenvalues and the variation of τ gives one additional negative eigenvalue.

Starting from an initial guess, we carry out a stationary-point search and thereby simultaneously optimize the bead positions, {x(0), …, x(N−1)}, and τ. In the normal regime,52 we use a standard single-index saddle point optimization algorithm.100 However, due to the nature of the instanton in the inverted regime as a higher-index saddle point, we have to adapt the optimization scheme slightly. Finding such higher-index saddle points can sometimes be a very demanding task. For instance, standard root-finding algorithms such as MINPACK’s hybrd and hybrj routines (based on a modified Powell method) as implemented, for example, in scipy can be used. Note that these approaches locate stationary points of any index and thus may not find the instanton unless a very good initial guess is given. However, this is made simpler as we exactly know the index of the saddle point we are seeking. This enables us to use eigenvector-following methods modified so as to reverse the signs of forces corresponding to the first N1D + 1 modes.101,102 One can alternatively modify a single-index saddle-point optimizer such as that introduced in Ref. 100 by reversing all but the first of these modes. This is the approach used to optimize the instantons in Sec. V.

One might worry that the number of high-index stationary points is often larger than the number of minima and single-index saddle points as was seen in Ref. 102 for studies of Lennard–Jones clusters. However, in our case, the stationary points of the action have the direct physical interpretation of two classical trajectories, which join together into a periodic orbit. At least for the simple systems we have studied, it is clear that there is only one periodic orbit that can exist and thus only one stationary point of the action. In fact, we ran root-finding algorithms (which optimize to any stationary point, regardless of its index) starting from randomly chosen initial conditions and did not find any other stationary points of the action. We, therefore, conclude that there is no particular problem in locating ring-polymer instantons in the inverted regime.

In practice, we make a sophisticated initial guess of the instanton geometry using our knowledge that the hopping point is located at the crossing seam and about the general shape of the instanton in the inverted regime. In addition, we can initiate the optimization at a high temperature, where the instanton is known to be collapsed at the minimum-energy crossing point, and then systematically cool down the system, which ensures an excellent initial guess for each step. In this way, the instanton optimization in the inverted regime can be made just as numerically efficient and reliable as in the normal regime.

Once the instanton orbit is found, the derivatives of the action with respect to x′, x″, and τ, which are required for the prefactor, can be evaluated in terms of the bead positions, potentials, gradients, and Hessians using the formulas given in the Appendix of Ref. 52. This allows the rate to be computed directly using the formula in Eq. (11).

The analytic-continuation procedure ensures that we retain the correct behavior in the high-temperature classical limit and also that instanton theory continues to give the exact result for a system of two linear potentials.51 Here, we shall test the numerical implementation for a multidimensional and an anharmonic system to check that it remains well behaved. We chose to apply our method to the same two model systems as Lawrence and Manolopoulos in their work on the extrapolation of Wolynes theory into the inverted regime.62 

The first model system under consideration is the spin-boson model at T = 300 K defined by the potentials
V0(x)=j=1D12mωj2(xj+ζj)2,
(28a)
V1(x)=j=1D12mωj2(xjζj)2ε,
(28b)
where ζj=cj/mωj2 and
ωj=ωctan(j12)π2D,  cj=mΛ2Dωj,
with reorganization energy Λ = 50 kcal mol−1 and characteristic frequency ωc = 500 cm−1. This has a discretized spectral density in D degrees of freedom,
J(ω)=π2j=1Dcj2mωjδ(ωωj),
(29)
which reproduces a Debye spectral density in the D → ∞ limit.103,104
The exact quantum golden-rule rate for this system with constant diabatic coupling, Δ, can be calculated by numerical integration of the flux correlation function,88 
kQM=Δ22iτiτeϕ(t)/dt,
(30)
with
ϕ(t)=iεt+4πJ(ω)ω21cosωttanh12βω+isinωtdω,
(31)
where the rate is independent of τ which can therefore be chosen in order to get a faster convergence of the time-integral.
A commonly used approach105 is to perform a stationary-phase approximation to the rate expression in Eq. (30) about the point t = −iτ such that ϕ′(−iτ) = 0 to give
kSP=Δ222πϕ(iτ)eϕ(iτ)/,
(32)
where the primes denote derivatives with respect to t.
In the case of the spin-boson model, the closed-form expressions for the actions of classical trajectories on each surface are known50 and the instanton rate can be calculated analytically.51,55 The derivation, starting from Eq. (14), in the inverted regime is completely analogous. Note that the semiclassical partition function is exact for this harmonic system, i.e., Z0=j=1D[2sinhβωj/2]1. The action at the stationary point in the spatial coordinates is given by
S(τ)=ετ+j=1D2mωjζj21coshωjτtanh12βωj+sinhωjτ,
(33)
which holds for both positive and negative τ. In the case of the spin-boson model, the prefactor in Eq. (14) can be shown to cancel with the reactant partition function.51 Therefore, the rate expression reduces to
kSCI=2πΔ22d2Sdτ212eS(τ)/,
(34)
which should be evaluated at a value of τ found numerically to be the stationary point of the action [Eq. (33)]. This expression, however, coincides exactly with the stationary-phase approximation given in Eq. (32), as we identify S(τ) ≡ ϕ(−iτ). This therefore shows that, for the spin-boson model, the analytically continued instanton theory is equivalent to the stationary-phase approximation in both the normal and inverted regimes.
In order to demonstrate that the ring-polymer instanton optimization and rate calculation are numerically stable, we carried out calculations using Eq. (11) and a general multidimensional algorithm, which did not take the fact that we could solve the problem analytically into account. The bath was discretized into D = 12 degrees of freedom as in Ref. 62. We present all results as functions of the driving force, ε, and compare the computed rate constants with those of Marcus theory,
kMT(ε)=Δ2πβΛeβ(Λε)2/4Λ.
(35)
Note that Marcus theory is equal to classical TST for this system, to which instanton theory tends in the classical limit.51 The inverted regime is defined by ε/Λ > 1.

The semiclassical instanton rates are depicted in Fig. 6 where they can be compared with the exact result and with Marcus theory. As expected, instanton theory gives very accurate rates, at least for this system, as we have explained that it matches the results obtained by the stationary-phase expression, which, in turn, is known to be in excellent agreement with the exact rate for this model.106 Furthermore, Fig. 6 confirms that nuclear quantum effects in the inverted regime can be much larger than those in the normal regime,44 causing a dramatic orders-of-magnitude increase in the rate compared with the classical Marcus theory. It is for this reason that we considered it of particular importance to develop the practical instanton approach described in this paper.

FIG. 6.

Marcus theory, semiclassical instanton (SCI), and exact quantum rates for a twelve-dimensional spin boson model. Results are presented as a function of the driving force relative to the classical Marcus theory rate of the symmetric (ε = 0) system.

FIG. 6.

Marcus theory, semiclassical instanton (SCI), and exact quantum rates for a twelve-dimensional spin boson model. Results are presented as a function of the driving force relative to the classical Marcus theory rate of the symmetric (ε = 0) system.

Close modal

Table I shows how the results converge with the total number of ring-polymer beads, N, for systems in the normal, activationless, and inverted regimes. It can be seen that when the beads are split among the two potentials according to the “natural” ratio N1/N0 = |τ1|/τ0, the rate converges very quickly. Even with only N = 128 beads, all rates are found to be converged to within 2.5% of the stationary-phase result. However, in general, τ is not known prior to the optimization. Hence, we also show the rates optimized with an equal split of the beads, N0 = N1. Although the instanton rates approach the stationary-phase results slower compared to the rates calculated with a natural ratio, convergence is again reached in all cases. Consequently, a proficient initial guess makes convergence faster but is not required to obtain accurate results. Furthermore, a simple approach suggests itself in which one could use the optimized τ value obtained from an instanton search with a low number of beads. The split of the beads can then be adjusted accordingly for more accurate calculations performed with a larger number of beads.

TABLE I.

Numerical results for the reaction rates of the spin-boson model (parameters defined in Sec. V A) computed using various methods given relative to the Marcus rate for the same system as log10[k(ε)/kMT(ε)]. The values of τ given are determined from the calculation of the stationary-phase expression. We optimize the instanton using two approaches for splitting the N beads into two sets, one with the natural bead ratio (to the nearest integers) defined by N1/N0 = rnat = |τ1|/τ0 (using τ obtained from the stationary-phase approximation) and the other with an equal split N1/N0 = r1/2 = 1. A cell with “−” indicates failure to find a stationary point with the correct index. In the limit N → ∞, the result tends to the stationary-phase (SP) approximation [Eq. (32)]. Exact rates are calculated by numerically integrating Eq. (30).

ε0.00.51.01.52.0
τ/βℏ0.50000.15890.0000−0.0430−0.0612
Nrnatrnatr1/2rnatr1/2rnatr1/2rnatr1/2
32 2.24 1.12 1.50 −0.25 0.60 1.37 − 7.00 − 
64 2.26 1.12 1.26 −0.26 0.04 1.31 2.17 7.01 − 
128 2.26 1.13 1.16 −0.26 −0.18 1.31 1.47 7.04 7.19 
256 2.26 1.13 1.14 −0.26 −0.24 1.32 1.36 7.05 7.08 
512 2.26 1.13 1.13 −0.26 −0.26 1.32 1.33 7.05 7.06 
1024 2.26 1.13 1.13 −0.26 −0.26 1.32 1.33 7.05 7.06 
2048 2.26 1.13 1.13 −0.26 −0.26 1.32 1.32 7.05 7.05 
SP 2.26 1.13 −0.26 1.32 7.05 
Exact 2.26 1.13 −0.25 1.31 7.05 
ε0.00.51.01.52.0
τ/βℏ0.50000.15890.0000−0.0430−0.0612
Nrnatrnatr1/2rnatr1/2rnatr1/2rnatr1/2
32 2.24 1.12 1.50 −0.25 0.60 1.37 − 7.00 − 
64 2.26 1.12 1.26 −0.26 0.04 1.31 2.17 7.01 − 
128 2.26 1.13 1.16 −0.26 −0.18 1.31 1.47 7.04 7.19 
256 2.26 1.13 1.14 −0.26 −0.24 1.32 1.36 7.05 7.08 
512 2.26 1.13 1.13 −0.26 −0.26 1.32 1.33 7.05 7.06 
1024 2.26 1.13 1.13 −0.26 −0.26 1.32 1.33 7.05 7.06 
2048 2.26 1.13 1.13 −0.26 −0.26 1.32 1.32 7.05 7.05 
SP 2.26 1.13 −0.26 1.32 7.05 
Exact 2.26 1.13 −0.25 1.31 7.05 

These results confirm that golden-rule instanton theory is not only as accurate but also as efficient in the inverted regime as it is in the normal regime. In fact, convergence to an almost perfect quantitative agreement is achieved even deep in the inverted regime, where the quantum rate has a 107-fold speed up due to tunneling.

In this section, we show that the approach is not restricted to harmonic systems but works just as well for anharmonic potentials, which is, of course, the main advantage of the instanton approach. We consider the predissociation model previously studied in Refs. 48 and 62, which is not only anharmonic and asymmetric but also contains an unbound state. The potentials are given as
V0(x)=12mω2x2,
(36a)
V1(x)=Dee2α(xζ)ε,
(36b)
with reduced units m = 1, ω = 1, De = 2, α = 0.2, ζ = 5, β = 3, and = 1, whereas ε is varied. Both states are depicted in Fig. 7 for one particular choice of the driving force, ε. We present results for a range of values of driving force relative to the reorganization energy given by Λ = Dee2αζ, which is the key parameter for determining the crossover between normal and inverted regimes. The exact quantum golden-rule rate expressed in the eigenbases of reactant and product states, {ψ0λ,ψ1ν},48 initialized by a thermal distribution of reactant states is calculated using Eq. (1). Just as in the spin-boson example, we give all rates relative to the corresponding classical rate. For the predissociation model, with Δ taken to be constant, the classical limit is given by the one-dimensional classical golden-rule transition-state theory (cl-TST) rate [Eq. (18)], where Z0cl=1/βω.
FIG. 7.

The diabatic potential-energy curves for the one-dimensional predissociation model for a particular choice of driving force, ε, in the inverted regime. The reorganization energy, Λ, is also indicated.

FIG. 7.

The diabatic potential-energy curves for the one-dimensional predissociation model for a particular choice of driving force, ε, in the inverted regime. The reorganization energy, Λ, is also indicated.

Close modal

The results are depicted in Fig. 8, again showing a remarkable agreement between the exact and instanton rates with a maximal relative error of 0.1%. This confirms that even in the inverted regime, the semiclassical instanton method based on localized trajectories, as shown in Fig. 3(b), gives an excellent approximation to the standard golden-rule method, which is typically defined in terms of delocalized wave-function overlaps, as illustrated in Fig. 1(b), although they might look very different at first glance. The order-of-magnitude deviation of the classical rate for large ε emphasizes the remarkable relevance of nuclear tunneling effects in these systems, especially in the inverted regime, which can be almost perfectly captured with our semiclassical instanton approach.

FIG. 8.

Semiclassical instanton (SCI), exact quantum (QM), and classical TST rates for the predissociation model. Results from the various methods are presented as a function of the driving force, ε, relative to the classical rate of the ε = 0 system.

FIG. 8.

Semiclassical instanton (SCI), exact quantum (QM), and classical TST rates for the predissociation model. Results from the various methods are presented as a function of the driving force, ε, relative to the classical rate of the ε = 0 system.

Close modal

Similar in spirit to Ref. 62, analytical continuation of semiclassical instanton theory can be used to extend the method to the inverted regime leading to excellent results for the studied model systems from both approaches. Note, however, that although Lawrence and Manolopoulos were able to achieve similarly accurate results with their approach,62 it was necessary for them to design a special functional form to treat this dissociative system, whereas we could apply an identical algorithm to the case of the spin-boson model.

Besides the calculation of electron-transfer rates, we want to stress another interesting possible application of the method for which the predissociation model provides a simple example. Instead of artificially shifting two potential-energy surfaces in order to simulate different regimes of electron transfer, the shift between the two surfaces could be caused by a variable external field. For instance, we consider a system with a ground electronic state |0⟩, which is uncoupled to an excited electronic state |1⟩, and these potential-energy surfaces may be well separated in energy. We can then study the situation of the interaction of the uncoupled system with a light field with continuous-wave frequency ωex and interpret the golden-rule result as a photodissociation spectrum. The two electronic states will now be coupled by the electric dipole operator μ(x^) instead of the nonadiabatic coupling, Δ(x^). Hence, the golden-rule limit is equivalent to a linear response treatment of the weak-field interaction.

The reason why we can use the instanton method for this problem is because, like the rate, it is also described by Fermi’s golden rule. The simple connection between the definitions of the rate defined by Eq. (1) and the total photodissociation cross section in the weak-coupling limit initialized by a thermal equilibrium distribution becomes apparent when looking at the formula for the total photodissociation cross section starting from a thermal equilibrium distribution28,107,108
σtot(ωex)=πωexε0cλeβE0λZ0|μλν|2δ(E0λ+ωexE1ν)dE1ν,
(37)
where c is the speed of light, ε0 is the vacuum permittivity, and μλν=ψ0λ(x)*μ(x)ψ1ν(x)dx. Note that in our example of a scattering excited state, the energies E1ν are continuous. Therefore, we have replaced the sum in Eq. (1) by an integral and used energy-normalized continuum wave functions ψ1ν. Here, we shall assume the transition dipole moment to be constant, also known as the Condon approximation. Hence, the total cross section is directly related to golden-rule rate theory by
σtot(ωex)=ωex2ε0ckQM(ωex),
(38)
where the rate constant, kQM(ℏωex), is computed with the reactant potential shifted by the photon energy, i.e., V0(x) → V0(x) + ℏωex, and Δ is replaced by μ. The rate thus depends on the photon frequency in the sense that it shifts the potential-energy surfaces relative to each other. Similar expressions can be given for spontaneous or stimulated emission.107 

Replacing the exact rate in Eq. (38) by the instanton or the classical TST rate, we obtain the various approximate simulated spectra shown in Fig. 9. In this case, we used the predissociation model [Eqs. (36)] with a fixed value of ε = −2Λ, which then describes the typical case of an excited electronic-state potential V1 high above the ground-state potential V0. The deviation of the classical cross sections again illustrates the sizeable effect of quantum nuclear effects, this time on the line shape of optical spectra. On the other hand, semiclassical instanton theory reaches graphical accuracy with the exact result.

FIG. 9.

Total photodissociation cross sections for the predissociation model obtained by semiclassical instanton (SCI), exact quantum (QM), and classical TST calculations. The dissociative excited-state potential, V1, is given an asymptotic energy of −ε = 2Λ and is coupled to the ground-state potential, V0, by a continuous-wave weak field with frequency ωex.

FIG. 9.

Total photodissociation cross sections for the predissociation model obtained by semiclassical instanton (SCI), exact quantum (QM), and classical TST calculations. The dissociative excited-state potential, V1, is given an asymptotic energy of −ε = 2Λ and is coupled to the ground-state potential, V0, by a continuous-wave weak field with frequency ωex.

Close modal

In order to showcase what our method can contribute to the description of spectra, it is worth making a comparison with standard approaches used in quantum chemistry. The simplest and probably the most common method is to calculate the vertical excitation energy from the ground state minimum, which corresponds to Λ − ε in our model (specifically, 3Λ according to our choice of parameters). This gives a single, discrete line with zero width in the spectrum at ℏωex = Λ − ε, and with this approach, one completely disregards the statistics (or dynamics) of both states. This method can be improved by assuming a classical Boltzmann distribution in the ground state, which can, for example, be sampled by molecular dynamics simulations. By calculating vertical excitation energies from different sampled configurations, the natural width of the absorption bands can be revealed. This corresponds to our classical TST calculations. The instanton method improves this description by including quantum nuclear effects for both the ground and excited states. For example, by taking the zero-point energy of the reactant potential into account, the instanton results correctly exhibit a red shift of the absorption maximum relative to the classical prediction. Also, due to tunneling in the normal regime, the absorption band exhibits an earlier onset, whereas the equivalent process in the inverted regime causes a slower decay of the band. The speed up of the transition rate induced by quantum tunneling, therefore, directly translates into adding longer tails to both sides of the absorption spectrum.

We do, however, ignore the real-time dynamics within the wells such that our method thus probes only the time scale of the fastest decay of the wave packet, which imprints itself in the broadest features of the spectrum. These features form an envelope of constraint on the spectrum that will not be changed by any other dynamics leading to vibrational and rotational fine structure.109 Therefore, although we shall not be able to describe vibrational progressions with this approach, we expect to predict the correct envelope of the spectrum. We note that this example of transition to a dissociative state is thus a particularly favorable case for us.

We have extended the semiclassical instanton method for calculating Fermi’s golden rule into the Marcus inverted regime. It can be applied to multidimensional and anharmonic systems and gives a good approximation to the exact result even when nuclear quantum effects are important. The theory reduces to classical golden-rule TST in the high-temperature limit and hence Marcus theory when the potentials are harmonic, is exact for a system of two linear potentials, and is identical to the stationary-phase approximation in the case of the spin-boson model.

The main difference between the normal and inverted regimes is the form of the instanton periodic orbit, although in both cases, it is defined by the stationary point of the total action formed by joining two imaginary-time trajectories together. In the normal regime, both trajectories are minimum-action paths, which travel in positive imaginary time on the reactant or product potential-energy surfaces. However, in the inverted regime, the product trajectory is a maximum-action path, which travels in negative imaginary time. In both regimes, the energy and momentum are conserved when hopping from one state to the other, which occurs at a point where the two diabatic potentials are degenerate.

In order to locate the inverted-regime instanton within the ring-polymer formalism, we search for a high-index saddle point of the total action. We show that by using the knowledge we have about the expected number of negative eigenvalues as well as the approximate shape and location of the two trajectories, the algorithm can be made just as efficient as in the normal regime. Therefore, this approach can be used to calculate reaction rates across the whole range of electron-transfer reactions or for simulating spectral line shapes.

In contrast to closed-form expressions for the rate,22,87,110 which effectively require a global harmonic approximation for the potential-energy surfaces, the instanton approach locates the instanton without making any assumptions and takes only a local harmonic approximation about the tunneling path. All one has to provide to the algorithm are functions returning the potentials, gradients, and Hessians on the two diabatic potential-energy surfaces for a given nuclear configuration. Additionally, it only requires this knowledge about a rather small region located around the crossing point. This is another reason for the computational efficiency of the method and makes it conceptually easily applicable to molecular systems, even in conjuncture with on-the-fly ab initio electronic-structure codes. For further enhancements to the efficiency, machine-learning approaches could be applied.77 

Apart from the excellent agreement with the exact quantum rates for the model systems studied in this work, one of the main advantages of the method from a qualitative perspective is that it provides direct mechanistic insight into the reaction under investigation. In this respect, our method appealingly stands out from alternative approaches, which effectively extrapolate the data collected in the normal regime.62 The instanton orbit can be interpreted as the “dominant tunneling pathway” and identifies which nuclear modes are involved in the reaction. In cases where there are competing reactions, the instanton approach will identify the dominant mechanism. Comparison to the classical (Marcus or golden-rule TST) rate allows an easy identification of the role of quantum nuclear effects (tunneling and zero-point energy), which are expected to be particularly important in the inverted regime. Kinetic isotope effects can also be easily predicted, which often provide the easiest connection to experimental data.111 

A limitation of all instanton methods is that the semiclassical approximation is not valid for fluxional environments such as encountered by reactions in solution.79 This is because many different stationary points exist in this case, which cannot be treated by a steepest-descent integration as the Gaussians introduced by this approximation would have a significant overlap.

Nonetheless, a good understanding of the instanton paths has helped in the development of a number of methods based on path-integral sampling, which are applicable to reactions in solution.69,70,112 We hope that the information obtained on the instanton in this work will help derive novel path-integral-based rate theories, which can describe the inverted regime more rigorously.

The method described in this work is, however, well suited to be applied to complex chemical reactions in the gas phase, in clusters, on surfaces, and in solids. A wide range of processes involving a transition between two electronic states can be studied in this way, so long as the coupling falls in the golden-rule regime. This encompasses not only electron-transfer reactions but also spectroscopy, intersystem crossing, and electrochemical processes. Showing the capability of the method in such applications will be an integral part of future work.

This work was financially supported by the Swiss National Science Foundation through SNSF Project No. 175696.

1.
S. C.
Althorpe
,
N.
Ananth
,
G.
Angulo
,
R. D.
Astumian
,
V.
Beniwal
,
J.
Blumberger
,
P. G.
Bolhuis
,
B.
Ensing
,
D. R.
Glowacki
,
S.
Habershon
,
S.
Hammes-Schiffer
,
T. J. H.
Hele
,
N.
Makri
,
D. E.
Manolopoulos
,
L. K.
McKemmish
,
T. F.
Miller
III
,
W. H.
Miller
,
A. J.
Mulholland
,
T.
Nekipelova
,
E.
Pollak
,
J. O.
Richardson
,
M.
Richter
,
P. R.
Chowdhury
,
D.
Shalashilin
, and
R.
Szabla
,
Faraday Discuss.
195
,
311
(
2016
).
2.
J. E.
Lawrence
and
D. E.
Manolopoulos
,
Faraday Discuss.
221
,
9
(
2020
).
3.
J. C.
Tully
,
J. Chem. Phys.
137
,
22A301
(
2012
).
4.
B. F. E.
Curchod
and
T. J.
Martínez
,
Chem. Rev.
118
,
3305
(
2018
).
5.
R.
Crespo-Otero
and
M.
Barbatti
,
Chem. Rev.
118
,
7026
(
2018
).
7.
G. A.
Worth
,
M. A.
Robb
, and
B.
Lasorne
,
Mol. Phys.
106
,
2077
(
2008
).
8.
N.
Makri
,
Int. J. Quantum Chem.
115
,
1209
(
2015
).
10.
Multidimensional Quantum Dynamics: MCTDH Theory and Applications
, edited by
H.-D.
Meyer
,
F.
Gatti
, and
G. A.
Worth
(
Wiley-VCH
,
Weinheim
,
2009
).
11.
H.
Wang
,
D. E.
Skinner
, and
M.
Thoss
,
J. Chem. Phys.
125
,
174502
(
2006
).
12.
G. W.
Richings
,
I.
Polyak
,
K. E.
Spinlove
,
G. A.
Worth
,
I.
Burghardt
, and
B.
Lasorne
,
Int. Rev. Phys. Chem.
34
,
269
(
2015
).
13.
H.-D.
Meyer
and
W. H.
Miller
,
J. Chem. Phys.
70
,
3214
(
1979
).
14.
G.
Stock
and
M.
Thoss
,
Phys. Rev. Lett.
78
,
578
(
1997
).
15.
G.
Stock
and
M.
Thoss
,
Adv. Chem. Phys.
131
,
243
(
2005
).
16.
J. O.
Richardson
and
M.
Thoss
,
J. Chem. Phys.
139
,
031102
(
2013
).
17.
M. A. C.
Saller
,
A.
Kelly
, and
J. O.
Richardson
,
J. Chem. Phys.
150
,
071101
(
2019
); e-print arXiv:1811.08830 [physics.chem-ph];
[PubMed]
M. A. C.
Saller
,
A.
Kelly
, and
J. O.
Richardson
,
Faraday Discuss.
221
,
150
(
2020
); e-print arXiv:1904.11847 [physics.chem-ph].
18.
J. E.
Runeson
and
J. O.
Richardson
,
J. Chem. Phys.
151
,
044119
(
2019
); e-print arXiv:1904.08293 [physics.chem-ph].
19.
J. C.
Tully
,
J. Chem. Phys.
93
,
1061
(
1990
).
20.
J. E.
Subotnik
,
A.
Jain
,
B.
Landry
,
A.
Petit
,
W.
Ouyang
, and
N.
Bellonzi
,
Annu. Rev. Phys. Chem.
67
,
387
(
2016
).
21.
S.
Mai
,
P.
Marquetand
, and
L.
González
,
Wiley Interdiscip. Rev.: Comput. Mol. Sci.
8
,
e1370
(
2018
).
22.
I.
Rips
and
E.
Pollak
,
J. Chem. Phys.
103
,
7912
(
1995
).
23.
P. A. M.
Dirac
,
Proc. R. Soc. A
114
,
243
(
1927
).
25.
R.
Zwanzig
,
Nonequilibrium Statistical Mechanics
(
Oxford University Press
,
New York
,
2001
).
26.
A.
Nitzan
,
Chemical Dynamics in Condensed Phases: Relaxation, Transfer, and Reactions in Condensed Molecular Systems
(
Oxford University Press
,
Oxford
,
2006
).
27.
E.
Fermi
,
Nuclear Physics
(
University of Chicago Press
,
Chicago
,
1974
).
28.
R.
Schinke
,
Photodissociation Dynamics
(
Cambridge University Press
,
Cambridge
,
1993
).
29.
Conical Intersections: Theory, Computation and Experiment
, edited by
W.
Domcke
,
D. R.
Yarkony
, and
H.
Köppel
(
World Scientific
,
Singapore
,
2011
).
30.
M. A.
Parker
,
Physics of Optoelectronics (Optical Science and Engineering)
(
CRC Press
,
Boca Raton
,
2005
).
31.
R. F.
Ribeiro
,
L. A.
Martínez-Martínez
,
M.
Du
,
J.
Campos-Gonzalez-Angulo
, and
J.
Yuen-Zhou
,
Chem. Sci.
9
,
6325
(
2018
).
32.
33.
R. A.
Marcus
,
J. Chem. Phys.
24
,
966
(
1956
).
34.
R. A.
Marcus
and
N.
Sutin
,
Biochim. Biophys. Acta
811
,
265
(
1985
).
35.
J. R.
Miller
,
L. T.
Calcaterra
, and
G. L.
Closs
,
J. Am. Chem. Soc.
106
,
3047
(
1984
).
37.
A.
Rosspeintner
,
B.
Lang
, and
E.
Vauthey
,
Annu. Rev. Phys. Chem.
64
,
247
(
2013
).
38.
T.
Koslowski
,
F.
Burggraf
,
S.
Krapf
,
T.
Steinbrecher
, and
C.
Wittekindt
,
Biochim. Biophys. Acta
1817
,
1955
(
2012
).
39.
K. M.
Pelzer
and
S. B.
Darling
,
Mol. Syst. Des. Eng.
1
,
10
(
2016
).
40.
D.
Antoniou
,
J.
Basner
,
S.
Núñez
, and
S. D.
Schwartz
,
Chem. Rev.
106
,
3170
(
2006
).
41.
G.
Feher
,
J. P.
Allen
,
M. Y.
Okamura
, and
D. C.
Rees
,
Nature
339
,
111
(
1989
).
42.
R. M.
Lynden-Bell
,
Electrochem. Commun.
9
,
1857
(
2007
).
43.
D.
Chandler
, in
Classical and Quantum Dynamics in Condensed Phase Simulations
, edited by
B. J.
Berne
,
G.
Ciccotti
, and
D. F.
Coker
(
World Scientific
,
Singapore
,
1998
), Chap. 2, pp.
25
49
.
44.
P.
Siders
and
R. A.
Marcus
,
J. Am. Chem. Soc.
103
,
748
(
1981
).
45.
E.
Åkesson
,
G. C.
Walker
, and
P. F.
Barbara
,
J. Chem. Phys.
95
,
4188
(
1991
);
E.
Åkesson
,
A. E.
Johnson
,
N. E.
Levinger
,
G. C.
Walker
,
T. P.
DuBruil
, and
P. F.
Barbara
,
J. Chem. Phys.
96
,
7859
(
1992
).
46.
J.
Jortner
and
M.
Bixon
,
J. Chem. Phys.
88
,
167
(
1988
).
47.
P.
Siders
and
R. A.
Marcus
,
J. Am. Chem. Soc.
103
,
741
(
1981
).
48.
J. O.
Richardson
and
M.
Thoss
,
J. Chem. Phys.
141
,
074106
(
2014
); e-print arXiv:1406.3144 [physics.chem-ph].
49.
J.
Ulstrup
,
Charge Transfer Processes in Condensed Media
(
Springer-Verlag
,
Berlin
,
1979
).
50.
R. P.
Feynman
and
A. R.
Hibbs
,
Quantum Mechanics and Path Integrals
(
McGraw-Hill
,
New York
,
1965
).
51.
J. O.
Richardson
,
R.
Bauer
, and
M.
Thoss
,
J. Chem. Phys.
143
,
134115
(
2015
); e-print arXiv:1508.04919 [physics.chem-ph].
52.
J. O.
Richardson
,
J. Chem. Phys.
143
,
134116
(
2015
); e-print arXiv:1508.05195 [physics.chem-ph].
53.
J.
Mattiat
and
J. O.
Richardson
,
J. Chem. Phys.
148
,
102311
(
2018
); e-print arXiv:1708.06702 [physics.chem-ph].
54.
J.
Cao
,
C.
Minichino
, and
G. A.
Voth
,
J. Chem. Phys.
103
,
1391
(
1995
).
55.
J.
Cao
and
G. A.
Voth
,
J. Chem. Phys.
106
,
1769
(
1997
).
56.
C. D.
Schwieters
and
G. A.
Voth
,
J. Chem. Phys.
108
,
1055
(
1998
).
57.
S.
Ranya
and
N.
Ananth
, “
Multistate ring polymer instantons in the nonadiabatic limit
,” e-print arXiv:1908.03772 [physics.chem-ph].
58.
P. G.
Wolynes
,
J. Chem. Phys.
86
,
1957
(
1987
).
59.
P. G.
Wolynes
,
J. Chem. Phys.
87
,
6559
(
1987
).
60.
C. D.
Schwieters
and
G. A.
Voth
,
J. Chem. Phys.
111
,
2869
(
1999
).
61.
J. E.
Lawrence
,
T.
Fletcher
,
L. P.
Lindoy
, and
D. E.
Manolopoulos
,
J. Chem. Phys.
151
,
114119
(
2019
).
62.
J. E.
Lawrence
and
D. E.
Manolopoulos
,
J. Chem. Phys.
148
,
102313
(
2018
).
63.
S.
Habershon
,
D. E.
Manolopoulos
,
T. E.
Markland
, and
T. F.
Miller
III
,
Annu. Rev. Phys. Chem.
64
,
387
(
2013
).
64.
A. R.
Menzeleev
,
N.
Ananth
, and
T. F.
Miller
III
,
J. Chem. Phys.
135
,
074106
(
2011
).
65.
P.
Shushkov
,
J. Chem. Phys.
138
,
224102
(
2013
).
66.
A. R.
Menzeleev
,
F.
Bell
, and
T. F.
Miller
III
,
J. Chem. Phys.
140
,
064103
(
2014
);
[PubMed]
J. S.
Kretchmer
and
T. F.
Miller
III
,
Faraday Discuss.
195
,
191
(
2016
).
[PubMed]
67.
J. R.
Duke
and
N.
Ananth
,
Faraday Discuss.
195
,
253
(
2016
).
68.
X.
Tao
,
P.
Shushkov
, and
T. F.
Miller
III
,
J. Phys. Chem. A
123
,
3013
(
2019
).
69.
M. J.
Thapa
,
W.
Fang
, and
J. O.
Richardson
,
J. Chem. Phys.
150
,
104107
(
2019
); e-print arXiv:1811.05874 [physics.chem-ph].
70.
W.
Fang
,
M. J.
Thapa
, and
J. O.
Richardson
,
J. Chem. Phys.
151
,
214101
(
2019
); e-print arXiv:1910.04020 [physics.chem-ph].
71.
Y.
Litman
,
J. O.
Richardson
,
T.
Kumagai
, and
M.
Rossi
,
J. Am. Chem. Soc.
141
,
2526
(
2019
); e-print arXiv:1810.05681 [physics.chem-ph].
72.
P.
Bajaj
,
J. O.
Richardson
, and
F.
Paesani
,
Nat. Chem.
11
,
367
(
2019
).
73.
V.
Ásgeirsson
,
A.
Arnaldsson
, and
H.
Jónsson
,
J. Chem. Phys.
148
,
102334
(
2018
).
74.
J. O.
Richardson
,
C.
Pérez
,
S.
Lobsiger
,
A. A.
Reid
,
B.
Temelso
,
G. C.
Shields
,
Z.
Kisiel
,
D. J.
Wales
,
B. H.
Pate
, and
S. C.
Althorpe
,
Science
351
,
1310
(
2016
).
75.
J. B.
Rommel
,
Y.
Liu
,
H.-J.
Werner
, and
J.
Kästner
,
J. Phys. Chem. B
116
,
13682
(
2012
).
76.
J. O.
Richardson
,
J. Chem. Phys.
148
,
200901
(
2018
).
77.
G.
Laude
,
D.
Calderini
,
D. P.
Tew
, and
J. O.
Richardson
,
Faraday Discuss.
212
,
237
(
2018
); e-print arXiv:1805.02589 [physics.chem-ph].
78.
J. O.
Richardson
,
J. Chem. Phys.
144
,
114106
(
2016
); e-print arXiv:1512.04292 [physics.chem-ph].
79.
J. O.
Richardson
,
Int. Rev. Phys. Chem.
37
,
171
(
2018
).
80.
W. H.
Miller
,
S. D.
Schwartz
, and
J. W.
Tromp
,
J. Chem. Phys.
79
,
4889
(
1983
).
81.
M. C.
Gutzwiller
,
Chaos in Classical and Quantum Mechanics
(
Springer-Verlag
,
New York
,
1990
).
82.
W. H.
Miller
,
J. Chem. Phys.
55
,
3146
(
1971
).
83.
C. M.
Bender
and
S. A.
Orszag
,
Advanced Mathematical Methods for Scientists and Engineers
(
McGraw-Hill
,
New York
,
1978
).
84.
M. J.
Ablowitz
and
A. S.
Fokas
,
Complex Variables: Introduction and Application
,
2nd ed.
, Cambridge Texts in Applied Mathematics (
Cambridge University Press
,
Cambridge
,
2003
).
85.

This is effectively the same result as was obtained in Ref. 55 from an analysis of the “barrier partition function,” although in that paper, the extra approximation that C0C1CZ02 was also made, which is exact for the spin-boson model but not in general.

86.
W. H.
Miller
,
J. Chem. Phys.
62
,
1899
(
1975
).
87.
H.
Nakamura
,
Nonadiabatic Transition: Concepts, Basic Theories and Applications
,
2nd ed.
(
World Scientific
,
Singapore
,
2012
).
88.
U.
Weiss
,
Quantum Dissipative Systems
,
4th ed.
(
World Scientific
,
Singapore
,
2012
).
89.

The argument can easily be extended for the archetypal one-dimensional crossed linear system defined by Vn(x) = κn(xx) for which Sn=mx2/2τnκn(x+x)τnκn2τn3/24m. Solving for the stationary point Sx+=0 gives τ = βℏκ0/(κ0κ1), which is in agreement with Eq. (17). The other two equations give x = 0 and x+ = x as expected.

90.

In order for the rate to exist, we assume that at least one of the reactant and product has quasicontinuous states. If the states are truly continuous, one should replace the sums by integrals.

91.

Note that there will also be recurrences at later times for this finite-dimensional bound system, which cannot be captured by the steepest-descent approximation.

92.
R. P.
Feynman
,
Elementary Particles and the Laws of Physics
(
Cambridge University Press
,
Cambridge
,
1986
), Chap. The reason for antiparticles, pp.
1
60
.
93.

Imagine watching a movie backwards of a clock falling out of a window. It still has momentum pointing down, but it travels along a path upwards from the ground to the window and registers a negative elapsed time for this journey.

94.
L. D.
Landau
and
E. M.
Lifshitz
,
Mechanics
(
Butterworth-Heinemann
,
Oxford
,
1976
).
95.
C.
Zhu
and
H.
Nakamura
,
J. Chem. Phys.
102
,
7448
(
1995
).
96.
S.
Andersson
,
G.
Nyman
,
A.
Arnaldsson
,
U.
Manthe
, and
H.
Jónsson
,
J. Phys. Chem. A
113
,
4468
(
2009
).
97.
J. O.
Richardson
and
S. C.
Althorpe
,
J. Chem. Phys.
131
,
214106
(
2009
).
98.
J. B.
Rommel
,
T. P. M.
Goumans
, and
J.
Kästner
,
J. Chem. Theory Comput.
7
,
690
(
2011
).
99.

In the activationless regime, the value of τ at the stationary point will be 0, and thus, η1 = 0, which leads to problems in evaluating the ring-polymer action. However, this causes no problems as in this case, we know that the stationary point has all beads collapsed at the hopping point, which is also the minimum of V0, so no ring-polymer optimization is necessary. Derivatives of the action can then be evaluated analytically.

100.
J.
Simons
,
P.
Jerrgensen
,
H.
Taylor
, and
J.
Orment
,
J. Phys. Chem.
87
,
2745
(
1983
).
101.
D. J.
Wales
,
Energy Landscapes
(
Cambridge University Press
,
Cambridge
,
2003
).
102.
J. P. K.
Doye
and
D. J.
Wales
,
J. Chem. Phys.
116
,
3777
(
2002
).
103.
H.
Wang
and
M.
Thoss
,
J. Phys. Chem. A
107
,
2126
(
2003
).
104.
T. C.
Berkelbach
,
D. R.
Reichman
, and
T. E.
Markland
,
J. Chem. Phys.
136
,
034113
(
2012
).
105.
J. S.
Bader
,
R. A.
Kuharski
, and
D.
Chandler
,
J. Chem. Phys.
93
,
230
(
1990
).
106.

The stationary-phase result deviates from the exact rate by 4% in the activationless regime, but the errors for the other cases tested are well below 1%.

107.
D. J.
Tannor
,
Introduction to Quantum Mechanics: A Time-Dependent Perspective
(
University Science Books
,
Sausalito, CA
,
2007
).
108.
D. E.
Manolopoulos
and
M. H.
Alexander
,
J. Chem. Phys.
97
,
2527
(
1992
).
109.
E. J.
Heller
,
The Semiclassical Way to Dynamics and Spectroscopy
(
Princeton University Press
,
Princeton
,
2018
).
110.
S.
Jang
and
M. D.
Newton
,
J. Phys. Chem. B
110
,
18996
(
2006
).
111.
K.
Karandashev
,
Z.-H.
Xu
,
M.
Meuwly
,
J.
Vaníček
, and
J. O.
Richardson
,
Struct. Dyn.
4
,
061501
(
2017
).
112.
C. L.
Vaillant
,
M. J.
Thapa
,
J.
Vaníček
, and
J. O.
Richardson
,
J. Chem. Phys.
151
,
144111
(
2019
); e-print arXiv:1908.03419 [physics.chem-ph].