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.

## I. INTRODUCTION

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 propagation^{9–12} to quantum-classical mapping approaches^{13–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}

_{λν}between these states is given by the famous golden-rule formula

^{23,24}generalized for an initial thermal distribution,

^{25,26}

^{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 1992^{32} 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 theory^{43,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 semiclassical^{51–58} and imaginary-time path-integral approaches^{59–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 dynamics^{63} 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 dynamics^{66–68} and golden-rule quantum transition-state theory^{69,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 method^{51–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}

## II. INSTANTON THEORY IN THE NORMAL REGIME

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}

*n*∈ {0, 1} is the electronic-state index and x = (

*x*

_{1}, …,

*x*

_{D}) are the Cartesian coordinates of

*D*nuclear degrees of freedom. These nuclei move on the potential-energy surface

*V*

_{n}(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 $\Delta (x^)$ to give the total Hamiltonian in the diabatic representation,

^{43}

^{,}

^{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\lambda =E1\nu $ 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.

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.

### A. Correlation function formalism

^{43,48,80}

*Ĥ*

_{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*(i*z*) is an analytic function of *z* = *t* − i*τ* and as such the integral $\u222b\u2212\u221e\u221ec(iz)\u2009dz$ is independent of the contour of integration, and hence, the rate is independent of the choice of *τ*, at least within its range of validity.

_{i}to the final position x

_{f}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.

### B. Semiclassical approximation

^{81}generalized for imaginary-time arguments,

^{79,82}

_{n}(

*u*), which travels from x

_{n}(0) = x

_{i}to x

_{n}(

*τ*

_{n}) = x

_{f}in imaginary time

*τ*

_{n}. This trajectory is found as the path that makes the Euclidean action,

*S*

_{n}, stationary, and the action is defined as

*t*, employing the method of steepest descent.

^{83}This leads to the following expression for the golden-rule instanton rate:

^{51}

^{,}

*τ*, which describes the stationary point of the action, defined by $\u2202S\u2202x\u2032=\u2202S\u2202x\u2033=\u2202S\u2202\tau =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 relations

^{84}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,

*Z*

_{0}, should also be evaluated within a semiclassical approximation.

^{79}

^{85}

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)=mx\u0307n(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 $\u2202S\u2202\tau =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 *V*_{0}(x^{‡}) = *V*_{1}(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}

### C. Classical limit

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.

*x*

_{−}=

*x*

_{f}−

*x*

_{i}and $x+=12(xi+xf)$. The stationary points of the total action [Eq. (12)] can be found by searching first for the solution to $\u2202S\u2202x\u2212=0$, which gives

*x*

_{−}= 0, and then for the solution to $\u2202S\u2202\tau =0$ evaluated at

*x*

_{−}= 0, which requires that the hopping point,

*x*

_{+}=

*x*

^{‡}, obeys

*V*

_{0}(

*x*

^{‡}) =

*V*

_{1}(

*x*

^{‡}). Finally, $\u2202S\u2202x+=0$ requires that

*V*

_{n}(

*x*

^{‡}) is the derivative of the potential

*V*

_{n}(

*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 $\epsilon =V0(xmin(0))\u2212V1(xmin(1))$ and the (product) reorganization energy as $\Lambda =V1(xmin(0))\u2212V1(xmin(1))$, where $xmin(n)$ is the minimum of *V*_{n}(*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 regime^{51} but parallel in the inverted regime.

From Eq. (17), it can be seen that in the normal regime, where ∇*V*_{0}(*x*^{‡}) and ∇*V*_{1}(*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 ∇*V*_{0}(*x*^{‡}) = 0, the situation changes to *τ* = 0. In this paper, we shall consider only one type of inverted regime, ∇*V*_{1}(*x*^{‡})/∇*V*_{0}(*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, ∇*V*_{0}(*x*^{‡})/∇*V*_{1}(*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.

*V*

^{‡}=

*V*

_{0}(

*x*

^{‡}) =

*V*

_{1}(

*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 $\u222be\u2212\beta V0\delta (V0\u2212V1)\u2009dx$, 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 ∇*V*_{0}(*x*^{‡}) and ∇*V*_{1}(*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.

## III. INSTANTON THEORY IN THE INVERTED REGIME

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.

### A. Analytic continuation

The main complication for the derivation of instanton theory in the inverted regime is that the expression for the imaginary-time propagator *K*_{1} 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.

*Ĥ*

_{0}, labeled by the quantum number

*λ*and the eigenstates, $\psi 1\nu $, of the product Hamiltonian,

*Ĥ*

_{1}, labeled by

*ν*.

^{90}The flux correlation function [Eq. (5)] expanded in the energy basis is thus 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 $\u2211\nu e\u2212\tau nEn\nu /\u210f$ 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.

_{max}, for all

*ν*and x. Therefore, the absolute value of any term in the sum is necessarily less than $\Psi max2e\u2212\tau nEn\nu /\u210f$. 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

*K*

_{1}(x

_{i}, x

_{f},

*τ*

_{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\u0303(\tau )$ and the wave-function expansion for *K*_{n} 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\u0303(\tau )=c(\tau )$ only for 0 < *τ* < *βℏ*. Because $c\u0303(iz)$ and *c*(i*z*) 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 *K*_{1} is no longer uniformly convergent and thus *c*(*τ*) diverges. Interestingly, we find that the correlation function $c\u0303(\tau )$ 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\u0303(\tau )$ converges in this example and is analytic even for systems in the inverted regime where *τ* is negative.

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\u0303(\tau )$, 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*(i*z*) for positive *τ*, which, because $c(iz)=c\u0303(iz)$, must also be a good approximation to $c\u0303(iz)$ in this regime. Because $c\u0303(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 form^{2} 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.

### B. Analysis of the inverted-regime instanton orbit

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.

*S*

_{n}(x

_{i}, x

_{f},

*τ*

_{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, x

_{n}(

*u*), starts at x

_{n}(0) = x

_{i}and travels to x

_{n}(

*τ*

_{n}) = x

_{f}in imaginary time

*τ*

_{n}. This trajectory has a conserved energy,

*E*

_{n}, because the Hamiltonian is time-independent. We can, therefore, add a zero under the integral

*u*) and equal to

*E*

_{n}. Inserting this definition into Eq. (21a) leads to

_{n}is an infinitesimal displacement vector pointing along the path in the direction from x

_{i}to x

_{f}. 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 x

_{i}to x

_{f}in positive imaginary time

*τ*

_{n}> 0, the momentum, p

_{n}(

*u*), will point along the path in the direction from x

_{i}to x

_{f}. However, for a trajectory traveling x

_{i}to x

_{f}in negative imaginary time

*τ*

_{n}< 0, the momentum, p

_{n}(

*u*), will point in the opposite direction, i.e., along the path in the direction from x

_{f}to x

_{i}.

From these equations, we can determine that *S*_{n}(x_{i}, x_{f}, *τ*_{n}) ≡−*S*_{n}(x_{i}, x_{f}, −*τ*_{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, p_{n}, points in the opposite direction from the particle flow, i.e., it is antiparallel to dx_{n}. In particular, if the zero potential is chosen below the instanton energy (for instance, at the reactant minimum), making *E*_{n} ≥ 0, then *S*_{n} will be positive when *τ*_{n} > 0 and negative when *τ*_{n} < 0. Therefore, whereas *S*_{0} remains positive just as in the normal regime, in the inverted regime, the value of *S*_{1} becomes negative.

^{81}

_{n}at the end marked x′ (and likewise for double primes), i.e., $p0\u2032=mx\u03070(0)$, $p0\u2033=mx\u03070(\tau 0)$, $p1\u2032=mx\u03071(\tau 1)$, and $p1\u2033=mx\u03071(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\u2032$ = $p1\u2032$, $p0\u2032\u2032$ = $p1\u2032\u2032$, and

*E*

_{0}=

*E*

_{1}. 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 p

_{0}and p

_{1}must be vectors with the same direction and the same magnitude. Because the third condition requires the energies of the trajectories to match,

*E*

_{0}=

*E*

_{1}, 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

*V*

_{0}(x

^{‡}) =

*V*

_{1}(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, x_{n}(*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(\tau 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.

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.

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, *C*_{0} and *C*_{1} 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 x_{1} is a maximum-action path, and thus, all eigenvalues of the matrix $\u2212\u22022S1\u2202x\u2033\u2202x\u2032$ are negative. *C*_{1}, 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 *C*_{1}, 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 x_{1} is a maximum-action path, it will have an even higher index in the space of paths as explained in Sec. IV.

### C. Hamilton–Jacobi formulation

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.

*S*

_{1}in a slightly different way to give

*s*, which defines a metric in the configuration space such that (d

*s*)

^{2}= dx · dx. This gives a line integral along the path, defined in such a way that ∫d

*s*is the path length. Note that

*W*

_{n}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.,

*W*

_{n}(x

_{i}, x

_{f},

*E*

_{n}) =

*W*

_{n}(x

_{f}, x

_{i},

*E*

_{n}). The relations

*S*

_{n}= ±

*W*

_{n}+

*E*

_{n}

*τ*

_{n}, where the sign is chosen to match the sign of

*τ*

_{n}, can thus be viewed as Legendre transforms obeying the conditions $\u2202Sn\u2202\tau n=En$ and $\u2202Wn\u2202En=\u2212|\tau n|$. It also follows that $\u2202Sn\u2202xi=\xb1\u2202Wn\u2202xi$ and $\u2202Sn\u2202xf=\xb1\u2202Wn\u2202xf$.

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

*W*≡

*W*(x′, x″,

*E*) appropriately in the normal and inverted regimes, the total action [Eq. (12)] of the instanton can also be written as

*E*=

*E*

_{0}=

*E*

_{1}is a function of (x′, x″,

*τ*) via $En=\u2202Sn\u2202\tau n$ or

*β*is a function of (x′, x″,

*E*) via $\beta \u210f=\u2212\u2202W\u2202E$.

*W*follows in each case from the required relationship in Eq. (26). Differentiating Eq. (26) using the chain rule, we find $\u2202S\u2202x\u2032=\u2202W\u2202x\u2032$ and $\u2202S\u2202x\u2033=\u2202W\u2202x\u2033$, 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 $\beta \u210f=\u2212\u2202W\u2202E$ 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, *W*_{n}, and decreases the other. In the normal regime, the derivative of the total abbreviated action, *W* = *W*_{0} + *W*_{1}, vanishes at the hopping point because here $\u2202W0\u2202x\u2032=\u2202W0\u2202x\u2033=\u2212\u2202W1\u2202x\u2032=\u2212\u2202W1\u2202x\u2033=p\u2021$, where *p*^{‡} = *p*_{0}(*x*^{‡}, *E*) = *p*_{1}(*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* = *W*_{0} − *W*_{1} and its derivative also vanishes at the hopping point because here $\u2202W0\u2202x\u2032=\u2202W0\u2202x\u2033=\u2202W1\u2202x\u2032=\u2202W1\u2202x\u2033=\u2212p\u2021$.

Furthermore, it is possible to show in this formulation that in the normal regime, $\u22022W\u2202x\u2032\u2202x\u2032=\u22022W\u2202x\u2033\u2202x\u2033=2m[\u2207V0(x\u2021)\u2212\u2207V1(x\u2021)]/p\u2021>0$, whereas in the inverted regime, $\u22022W\u2202x\u2032\u2202x\u2032=\u22022W\u2202x\u2033\u2202x\u2033=\u22122m[\u2207V0(x\u2021)\u2212\u2207V1(x\u2021)]/p\u2021<0$, and in both cases, $\u22022W\u2202x\u2032\u2202x\u2033=0$. Therefore, the instanton is defined by a minimum of *W*(*x*′, *x*″) in the normal regime, but a maximum in the inverted regime.

### D. Interpretation of the mechanism

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 e^{−S/ℏ} = e^{−βE} e^{−W/ℏ}, 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 *W*_{n} ≥ 0, this is clearly the case in the normal regime, where *W* = *W*_{0} + *W*_{1}. However, in the inverted regime, where *W* = *W*_{0} − *W*_{1}, we will need to show that *W*_{0} ≥ *W*_{1}. This can be justified, at least for a system similar to the schematic in Fig. 3(b), for which x_{1} follows the part of the same path as x_{0}, although it is expected to hold more generally. Using the fact that *V*_{0}(x) ≥ *V*_{1}(x) and thus *p*_{0}(x, *E*_{0}) ≥ *p*_{1}(x, *E*_{1}) 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 *W*_{n} in Eq. (25) that *W*_{0} ≥ *W*_{1} 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 be^{51} $kSCI/kcl=exp[\beta 3\u210f2\kappa 02\kappa 12/24m(\kappa 0\u2212\kappa 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 e^{−W/ℏ} 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\u2212(W0\u2212W1)/\u210f$, where the *W*_{n} 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 x_{1} trajectory follows the part of the same path as x_{0}.

## IV. RING-POLYMER INSTANTON FORMULATION

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 *N*_{n} equally spaced imaginary-time intervals of length *η*_{n} = *τ*_{n}/*N*_{n}. This same approach can be adapted for use in the inverted regime, where *τ*_{1} and, therefore, *η*_{1} will be negative, while *N*_{1}, of course, remains positive. The resulting *N* = *N*_{0} + *N*_{1} beads describe the ring polymer as shown in Fig. 5.

^{(0)}≡ x

^{(N)}in order to form a closed orbit. This is defined such that in the

*N*

_{0},

*N*

_{1}→ ∞ 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 *S*_{RP} 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 *N*_{0}:*N*_{1} 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\tau 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 *N*_{0}*D* positive and *N*_{1}*D* + 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),\u2026,x(N0\u22121)}$ and maximize it with respect to the beads ${x(N0+1),\u2026,x(N\u22121)}$ 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 *N*_{1}*D* + 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).

## V. APPLICATION TO MODEL SYSTEMS

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}

### A. Spin-boson model

*T*= 300 K defined by the potentials

^{−1}and characteristic frequency

*ω*

_{c}= 500 cm

^{−1}. This has a discretized spectral density in

*D*degrees of freedom,

*D*→ ∞ limit.

^{103,104}

^{88}

*τ*which can therefore be chosen in order to get a faster convergence of the time-integral.

^{105}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

*t*.

^{50}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=\u220fj=1D[2\u2061sinh\u2061\beta \u210f\omega j/2]\u22121$. The action at the stationary point in the spatial coordinates is given by

*τ*. 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

*τ*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.

*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,

^{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.

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 *N*_{1}/*N*_{0} = |*τ*_{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, *N*_{0} = *N*_{1}. 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.

ε/Λ
. | 0.0 . | 0.5 . | 1.0 . | 1.5 . | 2.0 . | ||||
---|---|---|---|---|---|---|---|---|---|

τ/βℏ
. | 0.5000 . | 0.1589 . | 0.0000 . | −0.0430 . | −0.0612 . | ||||

N
. | r_{nat}
. | r_{nat}
. | r_{1/2}
. | r_{nat}
. | r_{1/2}
. | r_{nat}
. | r_{1/2}
. | r_{nat}
. | r_{1/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.0 . | 0.5 . | 1.0 . | 1.5 . | 2.0 . | ||||
---|---|---|---|---|---|---|---|---|---|

τ/βℏ
. | 0.5000 . | 0.1589 . | 0.0000 . | −0.0430 . | −0.0612 . | ||||

N
. | r_{nat}
. | r_{nat}
. | r_{1/2}
. | r_{nat}
. | r_{1/2}
. | r_{nat}
. | r_{1/2}
. | r_{nat}
. | r_{1/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 10^{7}-fold speed up due to tunneling.

### B. Predissociation model

*m*= 1,

*ω*= 1,

*D*

_{e}= 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 Λ =

*D*

_{e}e

^{2αζ}, 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, ${\psi 0\lambda ,\psi 1\nu}$,

^{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/\beta \u210f\omega $.

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.

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 $\mu (x^)$ instead of the nonadiabatic coupling, $\Delta (x^)$. Hence, the golden-rule limit is equivalent to a linear response treatment of the weak-field interaction.

^{28,107,108}

*c*is the speed of light,

*ε*

_{0}is the vacuum permittivity, and $\mu \lambda \nu =\u222b\psi 0\lambda (x)*\mu (x)\psi 1\nu (x)\u2009dx$. Note that in our example of a scattering excited state, the energies $E1\nu $ are continuous. Therefore, we have replaced the sum in Eq. (1) by an integral and used energy-normalized continuum wave functions $\psi 1\nu $. 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

*k*

_{QM}(

*ℏω*

_{ex}), is computed with the reactant potential shifted by the photon energy, i.e.,

*V*

_{0}(x) →

*V*

_{0}(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 *V*_{1} high above the ground-state potential *V*_{0}. 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.

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.

## VI. CONCLUSIONS

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.

## ACKNOWLEDGMENTS

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

## REFERENCES

*Multidimensional Quantum Dynamics: MCTDH Theory and Applications*

*Nonequilibrium Statistical Mechanics*

*Chemical Dynamics in Condensed Phases: Relaxation, Transfer, and Reactions in Condensed Molecular Systems*

*Photodissociation Dynamics*

*Conical Intersections: Theory, Computation and Experiment*

*Physics of Optoelectronics (Optical Science and Engineering)*

*Classical and Quantum Dynamics in Condensed Phase Simulations*

*Charge Transfer Processes in Condensed Media*

*Quantum Mechanics and Path Integrals*

*Chaos in Classical and Quantum Mechanics*

*Advanced Mathematical Methods for Scientists and Engineers*

*Complex Variables: Introduction and Application*

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 $C0C1\u2248CZ02$ was also made, which is exact for the spin-boson model but not in general.

*Nonadiabatic Transition: Concepts, Basic Theories and Applications*

*Quantum Dissipative Systems*

The argument can easily be extended for the archetypal one-dimensional crossed linear system defined by *V*_{n}(*x*) = *κ*_{n}(*x* − *x*^{‡}) for which $Sn=mx\u22122/2\tau n\u2212\kappa n(x+\u2212x\u2021)\tau n\u2212\kappa n2\tau n3/24m$. Solving for the stationary point $\u2202S\u2202x+=0$ gives *τ* = *βℏκ*_{0}/(*κ*_{0} − *κ*_{1}), which is in agreement with Eq. (17). The other two equations give *x*_{−} = 0 and *x*_{+} = *x*^{‡} as expected.

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.

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.

*Elementary Particles and the Laws of Physics*

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.

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 *V*_{0}, so no ring-polymer optimization is necessary. Derivatives of the action can then be evaluated analytically.

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

*Introduction to Quantum Mechanics: A Time-Dependent Perspective*

*The Semiclassical Way to Dynamics and Spectroscopy*