Instanton formulation of Fermi’s golden rule in the Marcus inverted 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.³ 

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

Marcus was awarded with the Nobel Prize in chemistry in 1992³² 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.³⁵ 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,⁴³ 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.¹ 

In particular, it has been predicted that quantum tunneling effects can substantially speed up the rate in the inverted regime.⁴⁴ 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.⁴⁶ 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.⁴⁷ 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.⁵⁰ For instance, Fermi’s golden rule has been recast into a semiclassical instanton (SCI) formulation,⁵¹ 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,⁶² who analytically continued Wolynes’ rate expression,⁵⁹ 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⁶³ where the transferred electron is treated as an explicit particle,⁶⁴ as can be explained by a semiclassical analysis.⁶⁵ 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 

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

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.

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.

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.

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.

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 τ₀ = βℏ − τ and the other from x″ to x′ in imaginary time τ₁ = τ on the product state. This forms a closed path of total imaginary time τ₀ + τ₁ ≡ βℏ, known as the instanton.

Classical trajectories in imaginary time are described by ordinary classical mechanics but with an upside-down potential.⁸² 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.⁵¹ 

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 V₀(x^‡) = V₁(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.⁵¹ 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.⁵³ It has also been employed in quantifying the contributions from competing reaction mechanisms.⁷⁰ 

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.⁵¹ 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.

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.⁸⁷ 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 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⁵¹ but parallel in the inverted regime.

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

Equation (18) is, in fact, valid, not just in the normal regime but also for the inverted regime, that is, whether ∇V₀(x^‡) and ∇V₁(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.⁵¹ 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.⁵⁴ 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 τ₁ < 0 and τ₀ > βℏ in such a way that τ₀ + τ₁ ≡ βℏ 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 K₁ diverges for negative τ₁. 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.

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.

Inserting the definition for θ_λν into $c̃(τ)$ 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̃(τ)=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 K₁ 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.

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.⁹¹ 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.⁶² 

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.⁸⁴ 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 form² 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.

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₀ remains positive just as in the normal regime, in the inverted regime, the value of S₁ becomes negative.

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(τ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 τ₀ and then suddenly turns into a particle on the product state with the same energy and momentum where it travels for an imaginary time τ₁ 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.⁹² 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.⁹³ 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₀ and C₁ 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₁ is a maximum-action path, and thus, all eigenvalues of the matrix $−∂2S1∂x″∂x′$ are negative. C₁, 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₁, 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₁ 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.⁸¹ This leads to further insight into the behavior of the instanton in the inverted regime.

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,⁹⁴ and in the normal regime, both S₀ and S₁ are minima with respect to variation of the path. However, in the inverted regime, the product trajectory will be a maximum of S₁. Trajectories that bounce once have a conjugate point⁸¹ and so the associated W_n are single-index saddle points with respect to variation of the path.⁵¹ There is nothing different in the definition of W₁ in the inverted regime, so the product trajectory is also a single-index saddle point.

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₀ + W₁, vanishes at the hopping point because here $∂W0∂x′=∂W0∂x″=−∂W1∂x′=−∂W1∂x″=p‡$⁠, where p^‡ = p₀(x^‡, E) = p₁(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₀ − W₁ and its derivative also vanishes at the hopping point because here $∂W0∂x′=∂W0∂x″=∂W1∂x′=∂W1∂x″=−p‡$⁠.

Furthermore, it is possible to show in this formulation that in the normal regime, $∂2W∂x′∂x′=∂2W∂x″∂x″=2m[∇V0(x‡)−∇V1(x‡)]/p‡>0$⁠, whereas in the inverted regime, $∂2W∂x′∂x′=∂2W∂x″∂x″=−2m[∇V0(x‡)−∇V1(x‡)]/p‡<0$⁠, and in both cases, $∂2W∂x′∂x″=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 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₀ + W₁. However, in the inverted regime, where W = W₀ − W₁, we will need to show that W₀ ≥ W₁. This can be justified, at least for a system similar to the schematic in Fig. 3(b), for which x₁ follows the part of the same path as x₀, although it is expected to hold more generally. Using the fact that V₀(x) ≥ V₁(x) and thus p₀(x, E₀) ≥ p₁(x, E₁) 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₀ ≥ W₁ 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⁵¹  $kSCI/kcl=exp[β3ℏ2κ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 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−(W0−W1)/ℏ$⁠, 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₁ trajectory follows the part of the same path as x₀.

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 τ₁ and, therefore, η₁ will be negative, while N₁, of course, remains positive. The resulting N = N₀ + N₁ beads describe the ring polymer as shown in Fig. 5.

As in the normal regime, we define the instanton as a stationary point of S_RP with respect to the bead positions, {x⁽⁰⁾, …, x^(N−1)}, and τ simultaneously. Note that following our ring-polymer instanton approach for the normal regime,⁵² 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₀:N₁ 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 N₀D positive and N₁D + 1 negative eigenvalues. This can be understood by first considering the situation where the values of τ as well as x⁽⁰⁾ and $x(N0)$ are fixed at the hopping point. We would need to minimize the action with respect to the beads ${x(1),…,x(N0−1)}$ and maximize it with respect to the beads ${x(N0+1),…,x(N−1)}$ 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 η₁ < 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⁽⁰⁾, …, x^(N−1)}, and τ. In the normal regime,⁵² we use a standard single-index saddle point optimization algorithm.¹⁰⁰ 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₁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).

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.⁵¹ 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.⁶²