Nonadiabatic quantum transition-state theory in the golden-rule limit. I. Theory and application to model systems

We propose a new quantum transition-state theory for calculating Fermi's golden-rule rates in complex multidimensional systems. This method is able to account for the nuclear quantum effects of delocalization, zero-point energy and tunnelling in an electron-transfer reaction. It is related to instanton theory but can be computed by path-integral sampling and is thus applicable to treat molecular reactions in solution. A constraint functional based on energy conservation is introduced which ensures that the dominant paths contributing to the reaction rate are sampled. We prove that the theory gives exact results for a system of crossed linear potentials and also the correct classical limit for any system. In numerical tests, the new method is also seen to be accurate for anharmonic systems, and even gives good predictions for rates in the Marcus inverted regime.


I. INTRODUCTION
][10] This theory can be rigorously derived as an asymptotic approximation to the exact rate constant, 11,12 and the ringpolymer representation of the instanton can be employed to study polyatomic reactions in the gas phase. 4,13,147][28][29] The excellent accuracy achieved by RPMD rate theory can be explained by its connection to instanton theory, as the instanton can be shown to be equivalent to the dominant ring-polymer configuration sampled by the RPMD scheme. 4RPMD can of course also be used to simulate gas-phase reactions 30,31 where it has been shown that it gives a similar level of accuracy to that of the instanton. 32n this paper, we shall study reactions in the nonadiabatic limit, for which a typical case is electron transfer accompanied by a reorganization of the solvent environment. 33A simulation requires the consideration of at least two electronic states, and hence one cannot use the Born-Oppenheimer (BO) approximation to study these processes.Quantifying the rate of such reactions is challenging, but instrumental in understanding redox processes in chemical and biological systems. 34Nuclear quantum effects such as tunnelling can play a significant role in these reactions, [35][36][37] as has been experimentally confirmed by large room-temperature kinetic isotope effects. 38a) Electronic mail: jeremy.richardson@phys.chem.ethz.ch The celebrated Marcus theory [39][40][41] provides an approximate solution to the rate of an electron-transfer reaction by assuming that the nuclei can be described by classical statistical mechanics and that the resulting free-energy curves are harmonic.Although extensions have been proposed to extend Marcus theory beyond the classical approximation, many of them rely on an underlying assumption that the system can be described by a spin-boson model. 42Hence there is a need for a full-dimensional atomistic simulation approach which is applicable to study electron-transfer processes in liquid systems including nuclear quantum effects.We shall derive a new quantum transition-state theory to treat such processes, using a connection to semiclassical instanton theory to guide our derivation.
In this paper, we shall consider the rate of population transfer from the reactant electronic state, |0 , to the product electronic state, |1 , at a temperature defined by β = (k B T ) −1 .The total Hamiltonian of the system is given by Ĥ Each diabatic state, |n , has an associated nuclear Hamiltonian, Ĥn = p2 /2m + V n (x), which for the purpose of simplifying the presentation we shall first consider to be one-dimensional and extend the approach to multidimensional systems in Sec.V. We shall assume that the golden-rule limit applies, in which the coupling, ∆, between these electronic states is weak.Fermi's golden-rule rate constant is defined as where Z 0 = tr e −β Ĥ0 is the reactant partition function.
The transition-state partition function is defined as 44 = tr e − Ĥ0τ0/ δ( Ĥ0 − E) e − Ĥ1τ1/ δ( Ĥ1 − E) dE , arXiv:1811.05874v1 [physics.chem-ph]14 Nov 2018 where τ 0 ≡ τ and τ 1 ≡ β − τ are imaginary times corresponding to a particular splitting of the Boltzmann operator between the reactant and product states.We shall employ these definitions for τ n in terms of τ and β throughout this paper.Note that because this expression forces the reactant and product energies to be equal, the exact golden-rule rate is independent of the choice of τ . 45n previous work, we have rigorously derived a semiclassical instanton expression for the rate constant of this reaction. 12,46This approach includes the important quantum effects of tunnelling, delocalization and zeropoint energy within an → 0 asymptotic approximation.We have shown how the instanton can be written as a discretized ring polymer, 47 which makes numerical evaluation of the rate computationally efficient, and have successfully tested this method on an asymmetric system which generalizes the spin-boson model such that the free-energy curves become anharmonic. 48he instanton is formed of two imaginary-time classical trajectories, one on each state, joined together to form a closed ring.The length of time the trajectory travels on the reactant surface is τ 0 whereas its time on the product surface is τ 1 .The total time of the instanton path is thus β , as in other semiclassical theories which involve the Boltzmann distribution. 8,49In the instanton approach, the time τ = τ inst is chosen such that the energies of the two trajectories match.The instanton must therefore always be located in the vicinity of the crossing seam, where V 0 (x) = V 1 (x), and in fact the two trajectories join together smoothly on the crossing seam to form a periodic orbit.Instanton rate theory is exact for a system of crossed linear potentials, as defined in Appendix B. In the classical limit, it is well behaved and gives a steepestdescent approximation to the true classical rate. 46he main disadvantage of semiclassical instanton theories is that they cannot be directly applied to systems with explicit solvents.This is because there would be such a large number of instanton orbits, each corresponding to a different configuration of the solvent, that they could not all be located.Also, if the instantons are close to each other in configuration space, the individual steepest-descent approximations will overlap with each other and give a poor estimate of the configuration integral. 50Taking inspiration from quantum simulations of liquids, [51][52][53][54] one sees that a better method would sample over relevant path-integral configurations in order to obtain the rate.
Wolynes has developed a theory to predict the goldenrule rate from a path-integral Monte Carlo simulation, 55 which has been employed with atomistic simulations in a number of applications. 35,56The theory was derived as a steepest-descent approximation to the flux-flux correlation function with respect to the time variable only.This is very similar to the derivation of the so-called "quantum instanton" theory in the adiabatic limit described in Ref. 57.The same theory can also be obtained by analytic continuation of the imaginary free-energy. 58The resulting method can be computed in terms of an in-tegral over imaginary-time paths, x(t), with the cyclic boundary conditions, x(0) = x(β ).The paths are split over the two potential-energy surfaces according to the variable τ and have action S τ [x(t)], which is defined later in Eq. ( 9).The value of τ is chosen to minimize the integral.However, Wolynes theory is not directly connected to semiclassical instanton theory, and as well as sampling configurations near to the crossing seam, it also includes configurations far away.Like the adiabatic quantum instanton approach, 6 it does not tend to the correct classical limit in general anharmonic systems. 47We are thus interested in deriving a new path-integral sampling scheme that demonstrates a stronger connection to semiclassical instanton theory.
We would like to obtain a new approach which samples over paths subject to a constraint such that the integral is dominated by the instanton-like configurations near the crossing seam.Like other quantum transition-state theories, [1][2][3][4][5] it should not perform real-time propagation but be based only on a constrained imaginary-time pathintegral simulation.4][5] In Ref. 4, it was shown that the optimal dividing surface passes through the instanton and in this way the semiclassical instanton rate can be related to quantum transition-state theory and hence to RPMD.
In order to obtain such a method, we shall need to introduce a constraint functional which is obeyed by instanton configurations.Following Ref. 59, we suggest using a constraint which forces the energy of the reactant and product states to be equal.We therefore propose a new theory, called golden-rule quantum transition-state theory (GR-QTST), in which the rate constant is defined by Eq. ( 2) with the ansatz where, written in path-integral notation, 60 the effective action, φ(τ ), is defined as The dimensionless constraint functional, σ τ [x(t)], will be defined below such that it forces the energy of the reactant and product to match for every path in the ensemble.Our derivation will rely on the idea that the method should predominantly sample the instanton trajectory and fluctuations nearby.Therefore only paths in the region around the crossing seam, including the instanton configuration itself, should be located on the constraint hypersurface, σ τ [x(t)] = 0.In this paper, we shall obtain a definition for an appropriate constraint functional which has these required properties.Rather than variationally optimizing a constraint based on centroids and ring-polymer normal modes, we shall use the physically suggestive constraint that the energy of the reactants and products must be equal.We shall show how the resulting method tends to the correct classical limit for high temperatures or heavy masses.A ring-polymer discretization, which makes the method applicable to complex multidimensional systems, will also be presented.We shall investigate a number of model systems for which exact results are available for comparison and will apply the method to atomistic simulations in future work.

II. THEORY
In this section, we define the action functional, S τ [x(t)], and derive a constraint functional, σ τ [x(t)], with the properties required for the new GR-QTST method.In order to make the argument as clear as possible, we first treat a one-dimensional system.The advantage of using path-integral approaches is that they easily generalize to multidimensional systems and the necessary extra steps are outlined in Sec.V.
We shall employ the Lagrangian formalism of classical mechanics and its extension to quantum mechanics provided by Feynman's path-integral theory. 60The Lagrangian for imaginary-time dynamics in the electronic state |n , corresponding to the Hamiltonian Ĥn , is Note that the sign of the second term is opposite from that of the usual expression as it is known that imaginary-time dynamics is equivalent to classical dynamics with an inverted potential. 49he most important functional is the Euclidean action, which is defined as A path for which S n [x(t)] is stationary is called a trajectory and is denoted x(t). 61wo open-ended paths are combined into a ring, x(t), such that x(0) = x(β ), where imaginary time, t, runs from 0 to τ on the reactant state and from τ to β on the product.The total action is the sum of the individual actions of each of these paths, where the time integral is from 0 to τ on the reactant and from τ to β on the product.The coordinates where the path hops from one state to another are Solving a set of equations, = 0 yields a special trajectory, with x(t) = x inst (t) and τ = τ inst .This classical trajectory is called the instanton and it hops at where x ‡ is defined as the crossing point of the diabatic surfaces, such that V 0 (x ‡ ) = V 1 (x ‡ ).That the derivative of the action with respect to τ is zero ensures that the two trajectories of the instanton have the same energy. 46Our formulation of GR-QTST is motivated by a similar energy-matching constraint.
As will become apparent, it will be necessary to rewrite the Lagrangian in terms of a generalized coordinate, q ≡ q(x).The Euler-Lagrange equation which is satisfied by trajectories, q(t) ≡ q(x(t)), is for q = q and q = q, (11)   where the conjugate momentum of the new coordinate is and the right-hand side is Following standard manipulations discussed in many classical mechanics text books, it can be shown that the instantaneous energy, here defined as E n (q, q) = − ∂ Ln ∂ q q + L n , is conserved along a trajectory, q(t), as its time-derivative is zero if the Euler-Lagrange equation is obeyed.There is, however, no unique definition for a functional reporting on the energy of a non-classical path, q(t).One possible measure is the average of the instantaneous energy along the path: This is known as the thermodynamic estimator as it can also be derived from the derivative with respect to β of the partition function.There are however other functionals which can be designed to report on the energy of a path.The only requirement is that they return the correct energy for trajectories, but may give different results for a non-classical path.
At first, we attempted to use the thermodynamic estimator, Eq. ( 14), to define the constraint functional but found that although the instanton is located on the constraint hypersurface, it also allows non-instanton paths to dominate the sampling, as shown in Fig. 1.The thermodynamic energy contains a minus sign in front of the kinetic-energy term and thus this constraint may allow for sampling of highly-stretched paths which are unphysical.To avoid this problem, it is necessary to remove all terms involving ẋ.Thus we shall use classical-mechanical considerations to obtain a different energy functional in such a way that it does not change the result for the instanton trajectories, but which is not defined in terms of ẋ.][64] First we integrate Eq. (14b) by parts to get The first term in the integrand will be modified by applying the Euler-Lagrange equation, Eq. ( 11), even though it is not obeyed by a non-classical path.We shall also choose to define q(x) such that for the instanton trajectory q inst ≡ q(x inst ) = 0 and q inst ≡ q(x inst ) = 0, and thus the boundary terms can be neglected.This gives the definition of the virial energy functional, The two energy estimators are equal only for an instanton trajectory: i.e.Ēth n [x inst (t)] = Ēv n [x inst (t)].Their value is also equal to the instantaneous classical energy E n (q inst , qinst ) at all points along the instanton trajectory.
One can now see that it was necessary to change to the generalized coordinate system so that we could ignore the boundary term.In order to choose the form of q(x) correctly, we take into consideration all knowledge that we have about the instanton in general.We have already discussed that the energies of the two trajectories are equal, E 0 = E 1 , and are using this fact to define the constraint.We also know that the instanton hops between electronic states at the crossing point where V 0 (x) = V 1 (x), and this information can be used to obtain a definition for q(x).
One choice for the generalized coordinate might be q . This obviously goes to zero at the ends of the instanton as required.The definition also has a physical interpretation as this is often used as the reaction coordinate in classical simulations of electron-transfer reactions. 33,65,66However, in general, because ∂ 2 q ∂x 2 = 0 it will give a complicated form for ∂ Ln ∂q (see Eq. ( 13)) which involves the ẋ terms which we are trying to remove.A better choice is its linearized form: which is a first-order Taylor expansion around x + ≡ 1 2 (x + x ).This is also a suitable definition for the generalized coordinate as it clearly goes to 0 at the hopping points of the instanton as required.Note that this approach does not require knowledge of the location of the crossing seam which would not in general be known for a complex multidimensional system.This generalized coordinate transform gives a simple form for where For systems where V − (x) is linear, which includes the spin-boson model, s(x + ) will be defined at exactly the crossing point, x ‡ , for all paths.
We choose to define the constraint functional for GR-QTST in terms of the virial energy as The factor of β is included in the definition so as to make the functional dimensionless and the factor of 2 3 ensures that the method tends to certain important limits as shown in Appendices A and B.
The GR-QTST proposed in this paper gives the correct rate for all systems in the classical limit, as shown in Appendix A. It gives the exact quantum rate in the case of a system of crossed linear potentials.This is proved in Appendix B. This suggests that within the ansatz of our new method, we have chosen the correct functional form for σ τ [x(t)].Had we, for instance, allowed the boundary terms to remain or not performed the virial transformation, the proof would no longer hold.Note that our choice of coordinate transform is however not unique in general.There will be other definitions of q(x) with the correct properties which in the case of a linear model or classical limit reduce to give the correct result.In future work, we shall test other forms of the generalized coordinate to explore possible improvements to the method.

III. RING-POLYMER REPRESENTATION
In order to apply the GR-QTST method to complex systems for which analytical results are not available, one discretizes paths into N ring-polymer beads, x = {x 1 , . . ., x N }. 51Between neighbouring beads, there are N 0 intervals describing dynamics of the reactant system and N 1 ≡ N − N 0 of the product system.The functionals become functions of the positions of the ringpolymer beads, which we obtain by using finite differences for time derivatives and the trapezoidal rule on the integrals over time.In this representation, the total ac-tion, equivalent to Eq. ( 9), is given by with cyclic indices such that x 0 ≡ x N .With this definition, there are two simple discretization schemes.The first is defined to keep the N 0 : N 1 ratio fixed as τ is varied, which leads to unequal spring force constants around the ring polymer.This is the scheme used in our ring-polymer instanton approaches. 47,48The second is to vary the ratio N 0 : N 1 with τ so as to preserve the equality τ 0 /N 0 ≡ τ 1 /N 1 .In the limit of N → ∞, the same results will be obtained by either method although the convergence with respect to N appears to be slightly better for the first approach, especially when τ approaches its extremes of 0 or β .However, the second approach, which is similar to the discretization scheme proposed by Wolynes, 55 may be simpler to implement in a computationally efficient atomistic simulation and is the approach employed in this paper.
The ring-polymer version of the constraint function, Eq. ( 20), is 22) with the virial energy on each state defined by where There may of course be other finite-difference formulae which have better convergence with N → ∞, but we found that this one worked well for a reasonably small value of N .
The path integral in Eq. ( 6) becomes the following constrained configuration integral over the ring polymer: with the normalization constants Λ = 2πβ 2 /mN . 60he rate constant of the GR-QTST method is then defined by Eq. ( 2) with the approximation of Eq. ( 5) and the definitions Eqs. ( 22) and (24).The reactant partition function is defined in the usual way as with N 0 = N and N 1 = 0.The GR-QTST method is defined using the virial energy estimator.However, for comparison, we note that the thermodynamic estimator in ring-polymer form is This is used to define the thermodynamic constraint, which however should not be used in the GR-QTST ansatz as we shall show.Wolynes theory 55 is based on an unconstrained pathintegral simulation and defines the rate constant by Eq. ( 2) with the ansatz where τ * is chosen to maximize the effective action, φ u (τ * ).
A probability distribution of the sampled configurations is plotted in Fig. 1 for both the thermodynamic and virial constraints, as well as for an unconstrained version.As shown in the insets, even for values of τ far from τ inst , configurations obtained with the virial constraint resemble instanton-like configurations localized close to the crossing seam, x ‡ .However configurations sampled using either the thermodynamic constraint or without constraint can be found far away.As expected, the unconstrained simulation samples configurations which do not require that the energy of beads on the reactant and product surfaces match and which are thus clearly unphysical descriptions of the reaction.In each case, simulations carried out at τ = τ inst do sample the instanton, but the unconstrained simulation as well as simulations performed using the thermodynamic constraint also sample a much wider set of configurations which are not localized at the crossing seam.Therefore results using these two approaches, unlike when using the virial constraint, may be contaminated by unphysical configurations.
We next study the behaviour of the action with respect to τ .There exists a point, x , for each value of τ between 0 and β which minimizes the action, S τ (x ), while obeying the constraint.The constrained minimization is performed by introducing a Lagrange multiplier, µ, and solving the following (N + 1) equations: where ∇ denotes the derivative with respect to x.Note that for τ = τ inst , the solution x represents the ringpolymer instanton configuration, x inst (t), but is otherwise a non-classical path.V−(x+) at a given τ obtained from path-integral sampling with (a) no constraint, (b) the thermodynamic constraint, and (c) the virial constraint.(Note that here the probability is normalized for each value of τ such that pτ (E) dE = 1.)This system is the one-dimensional harmonic model defined in Sec.V B, with /λ = 0.5 and no bath modes.For this system, τinst/β ≈ 0.79.The regions which are rarely sampled are coloured yellow whereas the commonly sampled regions are blue.The crosses signify the instanton configurations.The insets show the path with the minimal action at the indicated values of τ and E. The x-axis of the insets is the nuclear configuration, x, and the y-axis shows the potential-energy surfaces and the instantaneous energy of the beads.
As an example, we choose a one-dimensional biased system in which both potential-energy surfaces are harmonic.In Fig. 2, we plot the value of the action, Eq. ( 21), after performing a constrained minimization for different values of τ .The function S τ (x ) is almost flat for this x FIG. 2. Minimized actions, Sτ (x ), with and without constraints for various values of τ are compared for the onedimensional harmonic system used in Figure 1.The minimizations were performed with N = 50.The black dot indicates the action corresponding to the ring-polymer instanton configuration computed using the same number of pathintegral beads.
system, and in fact, as shown in Appendix B, it becomes perfectly flat in the limiting case of the crossed linear model.When τ is equal to its value at the instanton, τ inst , the value of S τ (x ) is equal to the instanton action, S inst .
Had we minimized S τ (x) without the constraint, this function would be strongly dependent on τ , as shown by the dashed line.Even with the thermodynamic energy constraint the action still varies over orders of magnitude depending on τ .Both Figures 1 and 2 show clearly that only by using the virial energy constraint can we say that there is a strong connection to the instanton result for any value of τ .

IV. PATH-INTEGRAL MONTE CARLO IMPLEMENTATION
Here we employ a Monte Carlo importance sampling approach based on an extension of the ringpolymer transition-state sampling methods used in previous work. 4,27It would also be possible to compute the GR-QTST integral using thermodynamic integration and Metropolis Monte Carlo 68 or path-integral molecular dynamics schemes. 69This will be necessary to extend our method to the atomistic description of liquids, and we shall discuss such applications in future studies.

A. Importance sampling
After locating the constrained minimum, x , for a given τ value, as described in Sec.III, we shall perform a coordinate transformation similar to a normal-mode transformation.In order to do this, we start by expanding our total action as a Taylor series around this point: In what follows, we do not truncate the series, so no approximation is made.The transformation shall be defined by the orthogonal matrix U , such that the first column is a vector parallel to ∇σ τ (x ), i.e.U i1 = ∇ i σ τ (x )|∇σ τ (x )| −1 , and so all other column vectors are orthogonal to this and normalized. 70We shall make use of the fact that the elements of the vector U T ∇σ τ (x ) are all zero except the first, i.e.
We have not yet defined the other column vectors of U .They will be chosen so as to diagonalize a submatrix of A = U T ∇ 2 S τ (x )U such that A kk = δ kk a k for k = 2, . . ., N and k = 2, . . ., N .Because x is the constrained minimum, we know that a k > 0 for k = 2, . . ., N .The coordinate transformation is then defined as y ≡ U T (x − x ), which gives To obtain an efficient Monte Carlo importance sampling scheme for the integral in Eq. ( 24), we change to the new coordinates and then multiply and divide by a set of Gaussian distributions to give where here the action S τ (x) is the full form, Eq. ( 21).We have introduced the normalized distribution GR-QTST (SD) GR-QTST Plot showing the effective action for a one-dimensional harmonic system computed using GR-QTST, Eq. ( 32), and its steepest-descent (SD) approximation, Eq. ( 35).The results are scaled by the reactant partition function, Z0, to ease comparison with the multidimensional system-bath model presented in Sec.V B. The effective action, φinst, of the instanton method, Eq. ( 36), is marked with a black dot at τ = τinst.The system is the same as was used in Figures 1 and 2 and simulations were performed with N = 50.
and defined y 1 such that σ τ x(y 1 , y 2 , . . ., y N ) = 0. Finally, note that the term arising from the normalization of the delta function is given by In this form, the integral lends itself well to a Monte Carlo importance sampling evaluation with P (y) defining the sampling distribution.For each Monte Carlo step, y k is sampled from the normal distribution, y k ∼ N (0, /a k ), for k = 2, . . ., N .Then y 1 is set to y 1 which is found by a one-dimensional root search algorithm, starting from a guess which solves the first- For the simulations performed using the virial constraint, the guess obtained from solving this equation analytically was always very close to the correct root and thus there was no problem in obtaining y 1 numerically.In the N → ∞ limit and within statistical error, the Monte Carlo importance sampling scheme outlined above introduces no extra approximations into the GR-QTST method.It is the approach which we shall use to test the method on model systems for which analytical results are not available.However, one could also obtain a steepest-descent approximation to Eq. ( 32) by ignoring higher-order terms in Eqs. ( 29) and (30), which would be exact only for the linear system.This gives y 1 = 0 and makes the exponent independent of y k for k = 2, . . ., N .Therefore e −φSD(τ Using Eq. ( 35) in Eqs. ( 2) and ( 5) gives a steepest-descent approximation to the GR-QTST rate.
For comparison, under the ring-polymer instanton formalism, the effective action, φ inst , is defined at τ = τ inst as follows: where a k are the eigenvalues of A and S inst is the action corresponding to the instanton trajectory.The instanton rate is given by Eq. ( 2) with the ansatz Z ‡ ≈ β e −φinst/ .This formulation is clearly related to the steepest-descent version of GR-QTST, Eq. ( 35), and has exactly the same exponential factor for τ = τ inst .There is a minor difference in the prefactor, just as between a steepest-descent approximation to ring-polymer transition-state theory and standard ring-polymer instanton rate theory. 4e are now in a position to make numerical calculations of GR-QTST for a general system.We choose to illustrate the method using the example of the onedimensional harmonic system, for which results are shown in Fig. 3.
The effective actions of GR-QTST and its steepestdescent approximation are very similar.Although the curves of φ(τ ) and φ SD (τ ) have a different shape due to subtle differences, they are both seen to be almost independent of τ and within one decimal place give the same numerical value.It is relatively unimportant which value of τ is used to define the GR-QTST rate because φ(τ ) is approximately independent of this choice.However, for reasons explained in Appendix C, the ansatz of the new method is evaluated at τ * which maximizes φ(τ * ).
As expected, it is seen that the steepest-descent approximation, φ SD (τ ), is in excellent agreement with the instanton value, φ inst , for all values of τ .This is due to the similarity between Eqs. ( 35) and (36), and the fact that the exponent, S τ (x ), is exactly equal to S inst at τ = τ inst and almost constant with respect to τ as shown in Fig. 2.
Because of these two relationships, GR-QTST is therefore also strongly connected to instanton theory and predicts an effective action similar to φ inst at all values of τ .This is a positive result as instanton theory is rigorously derived 46 and gives accurate rate predictions for this type of system. 48In liquids, where many diabatic crossings exist, both the steepest-descent approximation to the GR-QTST formula and instanton theory will fail.However GR-QTST should still be applicable to liquids and is expected to sample all the crossing seams correctly.

B. One-dimensional anharmonic model
We have given an analytical proof that GR-QTST is exact for the crossed linear system in Appendix B. To test the new method on an anharmonic system, we study a one-dimensional model of the dissociation of a molecule via an electron transfer, which has been employed to benchmark rate calculations in previous work. 59,71The potential-energy surfaces are defined as The parameters are specified in reduced units such that = 1, mass m = 1 and frequency ω = 1.We study the system with parameters = 0, D = 2, α = 0.2 and x 0 = 5.This is a non-trivial test case because V − (x) has a nonlinear dependence on x such that s(x + ) appearing in Eq. ( 23) is not always equal to the crossing point, x ‡ .
For comparison with GR-QTST, we also present results from exact, classical, semiclassical and Wolynes theories.The exact rate is found using a wave-function expansion of Eq. ( 3) using the known eigenstates of the harmonic and Morse oscillators. 72,73Semiclassical instanton rates were computed using the formula presented in Ref. 47.The classical rate for this anharmonic system is given by Landau-Zener formula 74,75 in the golden-rule limit, 59,76 where x ‡ denotes the crossing point.This formula can be obtained from the classical limit of the instanton rate 46 and takes into account that reactive trajectories might have positive or negative values of momenta and will contribute equally. 76s shown in Figure 4, the GR-QTST method is accurate at a wide range of temperatures.In the high temperature limit, it tends rigorously to the exact result, as do classical Landau-Zener and the semiclassical instanton theories.For example, at β = 0.005, the GR-QTST method demonstrates very small error of about 0.2% in the rate constant.At lower temperatures, such as at β = 6, it still gives good predictions but shows slightly higher error of 8%.Wolynes theory is actually more accurate than GR-QTST for this system at lower temperatures.However, at higher temperatures, Wolynes theory fails to reproduce the classical limit, and demonstrates an error of over 41% at β = 0.005.In the classical limit, all beads should collapse to a point at the crossing between the two diabatic surfaces. 59However, because the integrand of Wolynes theory does not contain a constraint in the path integral, it samples many other configurations as well which contaminate the result 47 (see Fig. 1).Although this may be seen as a rather minor problem here, its unphysical behaviour in the classical limit could have much more drastic consequences in complex systems.FIG. 4. Rate constants calculated for the one-dimensional anharmonic system defined in Eqs.(37).GR-QTST rates obtained using Eq. ( 32) are compared against exact, Eq. ( 3), semiclassical instanton (SC), Eq. ( 36), classical Landau-Zener (LZ), Eq. ( 38), and Wolynes, Eqs. ( 27), theories.N = 100 was used for all path-integral methods.

V. MULTIDIMENSIONAL GENERALIZATION
In this section we generalize the theory presented above to treat a system with f nuclear degrees of freedom, x = (x 1 , . . ., x f ).Each electronic state is thus described by the Hamiltonian Ĥn = f j=1 p2 j /2m + V n (x).We have mass-scaled each degree of freedom such that they have the same effective mass, m.

A. Definition of the multidimensional theory
As before, we use generalized coordinates, q = {q 1 , . . ., q f }.The first one is defined similarly to the onedimensional case as where x + = 1 2 (x + x ) and x = x(0) and x = x(τ ).The other generalized coordinates are chosen to be orthogonal to q 1 (x) and defined as where v k form a set of vectors orthogonal to each other as well as to ∇V − (x + ).Using the Einstein summation rule, the Lagrangian for each electronic state can be written L n = 1 2 m ẋj ẋj + V n (x), and its derivatives are The Jacobian matrix for the transformation, J, has elements [J] kj = ∂ q k ∂xj and the inverse transformation is given by The thermodynamic energy functional for a path, x(t), is given by and by following the same steps as in Sec.II, the virial energy functional is found to be where ∂ Ln ∂q k q k = ∂ Vn ∂xj x j − s j (x + ) and s(x . This was derived using the Euler-Lagrange equation, Eq. ( 11), for each degree of freedom.As before, the boundary term can be shown to vanish for the instanton trajectory because the following relations are obeyed at the end points: q inst1 = q inst1 = 0 and q instk = q instk = 0 for k = 2, . . ., f .This follows from the proof in Ref. 46 that the instanton changes state on the crossing seam and at this point has velocity perpendicular to the seam, which is in the q 1 direction only.
All other formulae given in Sec.III and Sec.IV A are easily generalized by replacing the scalar x with the vector x and Λ by Λ f .If we also allow the coupling to become dependent on position, the prefactor ∆ 2 should be replaced by ∆(x )∆(x ), and moved inside the integral of Eq. ( 6). 59he GR-QTST rate tends correctly to the classical limit also for multidimensional systems.This is because for short paths, which occur in the classical limit, q 1 is simply V − (x + ).This gives σ τ (x) = βV − (x), which is the required constraint in this limit, 59 similarly to the onedimensional case outlined in Appendix A.

B. Application to system-bath model
The new GR-QTST method is applied to a multidimensional system-bath model, and the rates obtained are compared against existing methods.The potentialenergy surfaces of the system-bath model are defined as follows, where and V b describes the bath which couples x 1 to the other degrees of freedom: The bath is defined by the Ohmic spectral density, J(ω) = mγω e −ω/ωc , where mγ is the friction coefficient and ω c is the cut-off frequency.The bath is discretized into f − 1 modes as 77 where 27 The one-dimensional system (f = 1) has no bath and is thus frictionless, whereas all other systems are dissipative with the same damping, γ.Note that the Hamiltonian of this system-bath model is equivalent by unitary transformation to that of a spin-boson model with Debye spectral density. 78he system considered in our numerical tests has reorganization energy λ ≡ 2mξ 2 Ω 2 = 80 kcal mol −1 and the bias, , is varied from 0 to 2λ.The system frequency is Ω/2πc = 500 cm −1 and the bath cut-off frequency ω c /2πc = 500 cm −1 , where here c is the speed of light in vacuum.The friction coefficient is defined by γ = 0.001 E h , where E h is the Hartree energy.0][81] Rates are computed at a temperature of T = 300 K.Note that the rate for this system is independent of mass, m.
The classical rate constant, Eq. (A3), of the systembath model in the golden-rule limit is given by Marcus theory as 40 The exact Fermi's golden-rule rate was calculated using the analytical path-integral expression 77 for the flux-flux correlation function 45 and numerically integrated to the first plateau. 35n previous work, we have shown that instanton theory gives excellent predictions for the rate of the reaction in this model. 48However, the ring-polymer instanton method is not numerically applicable in the inverted regime because τ inst would fall outside of the range [0, β ].Nonetheless, one can use analytically obtained instantons to predict semiclassical rates in the inverted regime for a finite number of beads.To obtain an analytical expression for the instanton rates, we consider a coordinate transform of the multidimensional systembath Hamiltonian to a spin-boson Hamiltonian. 78,82Ω j are the normal-mode frequencies and ξ j are the displacements of the harmonic oscillators that result from this coordinate transformation.For N beads, the semiclassical instanton rate constant is given by the following formula: where inst is the action corresponding to an instanton trajectory discretized with N ring-polymer beads for τ = τ inst , and is defined as where Υ 83 In our calculations, we chose to distribute an equal number of beads on the reactant and product surfaces, i.e.N 0 = N 1 = N/2.To compute the rate in Eq. ( 52), one chooses τ inst numerically as the value of τ which solves d dτ S (N ) inst = 0. To calculate the GR-QTST rates, the Monte Carlo scheme we proposed in Sec.IV A is employed.The importance sampling scheme turns out to be extremely efficient for this problem and we are able to get an extremely small statistical error of less than 1% with about 10 4 random samples for this system-bath model with 8 degrees of freedom.
In Fig. 5, we show the form of the effective action, φ(τ ), to investigate its behaviour on adding extra degrees of freedom.In Appendix C we examined the effect of adding uncoupled bath modes to the linear crossing model.In this case where the bath is coupled to the system, the behaviour is similar in many ways.The effect on φ SD (τ ) is to simply shift it to slightly higher values but it remains approximately flat.However, adding extra degrees of freedom to the model has a larger effect on φ(τ ), which gets increasingly curved around its maximum with each additional degree of freedom.The value of τ * chosen to compute the GR-QTST rate is the one which maximizes these curves, and this value is found at approximately the same value as for the instanton, τ inst .At this point, the value of φ(τ * ) is very similar to that of the instanton, φ inst , and thus GR-QTST will predict a rate similar to that of instanton theory as we wanted.
Figure 6 shows the results for this multidimensional system-bath model for varying degrees of bias, .It is seen that the semiclassical and GR-QTST methods both The dependence of the effective action on τ for the system-bath model with /λ = 0.5 and in each case, a different number of dimensions, f .All calculations were performed using N = 50.GR-QTST results, computed using Eq. ( 32), are shown with statistical error bars.The result of the steepest-descent (SD) approximation, Eq. ( 35), is included only for f = 8 as for all multidimensional (f ≥ 2) cases, the SD approximations to φ(τ ) are approximately the same.The effective actions of the instanton method, φinst, are marked with dots at τ = τinst for each case, although these lie almost exactly on top of each other.All results are normalized by the partition function Z0 to ease comparison with the frictionless results shown in Fig. 3.
give excellent predictions for the rate both in the normal and inverted regimes, where > λ.Because in this case the reactant and product surfaces are harmonic, Wolynes theory gives the same result as semiclassical instanton theory.However, neither Wolynes theory nor ring-polymer instanton theory is numerically applicable to study nonadiabatic reactions in the inverted regime.It is therefore particularly interesting that GR-QTST gives accurate predictions even deep into this regime.
In the inverted regime, the optimal value of τ inst which makes the action stationary is greater than β .According to the ansatz of GR-QTST the rate should be evaluated at τ * = β , which gives the maximum φ(τ * ) which can be obtained.As shown in Figure 2, because S τ (x ) is quite flat with respect to τ and in fact exactly so in the linear case, the effective actions sampled will be similar to the action of the analytic instanton solution and thus the rate prediction, like that of instanton theory, will be a good approximation to the exact result.
We did not have to rely on analytic continuation of the function φ(τ ) outside the range [0, β ] in order to obtain this result because the GR-QTST rate is approximately independent of the value of τ .However, because it can be curved, especially with many degrees of freedom, it may be possible to obtain a better prediction using extrapolation to locate the optimal splitting τ * , as has been  51), semiclassical instanton (SC), Eq. ( 52), and GR-QTST.Wolynes theory is equal to SC for this case but is only numerically applicable in the normal regime, < λ.The reference, kMT(0), is the classical Marcus theory rate constant for = 0. SC rates were computed analytically using N → ∞, whereas the GR-QTST rates were obtained with N = 100 ring-polymer beads.
done previously from Wolynes theory in Ref. 71.As our function is flatter than the unconstrained form, the extrapolation will be more robust and less likely to give an unphysical prediction.
In order to analyse the accuracy of the various methods in more detail, the calculated rate constants are given in Table I.All results are reported relative to the classical Marcus theory rate, which is independent of f .
For the frictionless one-dimensional system, the exact rate is technically not defined as the quantum dynamics is coherent.Instanton theory and GR-QTST, like other transition-state theories, cannot describe this behaviour and thus for comparison, we give the value obtained by integrating the flux-flux correlation function up to the first plateau.We have checked that the friction of our system is strong enough to ensure that the revivals are destroyed and the exact rate is well defined for the multidimensional systems.About 8 degrees of freedom are needed to converge these rates to three significant figures.The exact calculations for this system indicate that the Fermi's golden-rule rates decrease slightly because of the coupling to the bath modes.This is a well-known effect in dissipative systems, in which friction hinders tunnelling. 77,80,84,85s has been observed in previous work, 48 the semiclassical instanton approximation has excellent accuracy in predicting the rates of the spin-boson model.This is due to the fact that as the potentials are harmonic, the   52), and GR-QTST calculations on the system-bath model defined in the text for different values of bias, , and number of degrees of freedom, f .Exact rate constants were computed by numerically integrating over the analytic expression for the flux-flux correlation function as far as the first plateau.In the table, we present the quantum factor k/kMT, where the Marcus theory rate constant was calculated using Eq. ( 51).The GR-QTST results are reported along with the Monte Carlo sampling error, which is given for the last digit in parentheses.
steepest-descent approximation is exact for all nuclear degrees of freedom.By comparison with semiclassical calculations with a finite number of beads, we also show that the discretization error is very small.In the normal regime, it is seen that N = 100 is enough to converge the results to three significant figures, whereas in the inverted regime, because of particularly large tunnelling effects, one needs N = 200.Rates obtained by GR-QTST are then compared with the instanton method for the same number of ringpolymer beads, N .It is seen that GR-QTST makes an error of less than 4% in each case in the normal regime, where < λ.However, deep in the inverted regime, where /λ = 2, the new method demonstrates slightly higher error of about 30%, but nonetheless continues to give a correct order-of-magnitude estimate of the rate even though the quantum effects are greater than a million.

VI. DISCUSSION
In this paper, we have described a new quantum transition-state theory for the calculation of nonadiabatic rates in the golden-rule limit.The current work only presents applications to low-dimensional model systems for which Monte Carlo importance sampling is simple and accurate to apply.However, in order to employ the new approach with atomistic simulations of electrontransfer reactions in solution, different sampling schemes will be required.Due to the quantum-classical correspondence of the path-integral formalism, many of the molecular dynamics sampling techniques developed for classical statistical mechanics can be used in an extended ring-polymer space to evaluate the imaginary-time path integration in a numerically efficient manner. 51,69,86For instance, the GR-QTST rate could be computed with biased path-integral molecular dynamics using a thermodynamic integration along the reaction coordinate σ τ (x).We shall explore these possibilities in future work.
In contrast to Wolynes theory, 55 the new method rigorously tends to the correct classical limit.For a system of crossed linear potential-energy surfaces uncoupled to a harmonic oscillator bath, it gives the exact Fermi's golden-rule rate.Furthermore, due to its connection to instanton theory, it also gives a good order-of-magnitude prediction for more general systems for a wide range of temperatures.It is particularly interesting that it is able to go beyond ring-polymer instanton theory and make accurate predictions for rates deep in the inverted regime.
Note that other forms of the ansatz, Eqs. ( 5), might also exist.
We originally considered Z ‡ ≈ β 0 e −φ(τ )/ dτ , which is similar to a rate formulation in terms of the Kubo transform. 87,88This was found to work well for one-dimensional systems, where φ(τ ) is almost flat, but was not as accurate for the case of multidimensional systems.In future work, we shall investigate variations of the constraint functional which could be used with an alternative ansatz.
2][3][4][5] This introduces the numerically challenging problem of locating the optimal dividing surface in each case which makes them difficult to apply in general to complex systems.One thus uses RPMD which is independent of this choice. 27,28,89In contrast to these approaches, the GR-QTST formulation employs a completely different constraint functional based on energy conservation which does not need to be variationally optimized.Because of this, we do not require dynamical methods like RPMD in this case.This work also points towards the development of a general nonadiabatic quantum transition-state theory which would be valid beyond the golden-rule approximation.In this more general regime, the exact rate constant can still be written as two Boltzmann factors separated by flux operators with an energy conservation requirement. 59The only difference is that hops between surfaces are also allowed within the Boltzmann operators.The imaginary-time path integrals for such problems are easily evaluated using matrix algebra. 90If an appropriate energy constraint functional can be devised, such a method would be able to describe a wide range of chemical reactions from those occurring in the golden-rule limit, where ∆ → 0, to those in the Born-Oppenheimer regime where ∆ is large.Some progress in this direction has been achieved already, 58,91-93 but these methods do not include energy constraints and are thus not directly comparable to our GR-QTST approach.The kinetically-constrained RPMD approach (KC-RPMD), [94][95][96] does employ an energy constraint to perform the path-integral sampling.The constraint is however applied only on the centroid which means that such a constraint surface will not pass through the instanton for asymmetric barriers and the instanton configuration will not be included in the pathintegral sampling.KC-RPMD is a dynamical method, which makes use of real-time trajectories to correct for recrossing effects, and is designed to predict rates for all coupling strengths, ∆.Dynamics are performed in an extended space including an auxiliary variable which reports on the formation of hops in the ring polymer.Due to all of these differences, it is hard to see exactly how KC-RPMD and GR-QTST are related.However, we can at least make some numerical comparisons.In Ref. 95, the authors report results for a symmetric system-bath model (called System A3). 97 For this system with intermediate friction, γ/Ω = 1, GR-QTST gives k/k cl = 0.048 (4), where k cl is the classical adiabatic TST rate. 80This is in good agreement with the exact result, 0.053, whereas the KC-RPMD result reported in that paper is about 5 times too large.For the same symmetric system, in the weaker dissipative regime, γ/Ω ≤ 0.1, where nuclear coherence effects become significant, 79,80 neither KC-RPMD nor GR-QTST is able to describe the enhancement of the exact rate.Instead, the GR-QTST rate becomes independent of fric-tion, which is in agreement with the result obtained by integrating the analytic flux-flux correlation function up to the first plateau.The KC-RPMD rate deep in the inverted regime ( = 0.236 E h ) for System B of Ref. 94  is also reported to be too large compared with the exact result by almost an order of magnitude.For the same system, the GR-QTST method gives k/k MT = 8.0(2) × 10 6 , demonstrating a smaller error compared to the exact result, 1.4 × 10 7 .
Other approaches extract the quantum rate constant from a real-time dynamical simulation of the nonadiabatic reaction.Examples include those based on extensions of RPMD to employ the mapping formalism, [98][99][100][101][102][103][104] surface hopping [105][106][107] or explicit electron dynamics in the position representation 108,109 as well a related instanton theory. 110There are also mixed quantumclassical dynamical methods 111,112 based on linearized and partially linearized path-integral simulations, [113][114][115] Bohmian dynamics 116 or the exact factorization of the time-dependent electron-nuclear wave function. 117An isomorphic Hamiltonian can be used to include quantum nuclear effects in nonadiabatic trajectory simulations. 118t is expected that in an accurate nonadiabatic RPMD simulation, the ring polymer will pass through transitionstate configurations similar to those sampled in the GR-QTST ensemble.We therefore expect that our study of GR-QTST will be of use in understanding and improving the dynamical methods discussed above.
As it is defined as the difference between energy functionals, the constraint functional, σ τ [x(t)] would be independent of an energy bias, V ‡ , and of shifting the x coordinates.This is due to having used the generalized coordinate q = q(x) to allow neglecting the boundary terms.
The constraint functional can then be written where Noting that changing the integration variables for the path integral to the Fourier coefficients γ (n) k introduces a Jacobian, J n , the integral can be written where Λ n = 2πτ n /mN n .The fluctuation terms are and Using the properties of a delta function and performing the integral over x + gives which we know by comparison with the exact pathintegral expression for a linear potential, 119 performing the remaining Gaussian integrals gives where In the limit of an infinite number of modes, and we find that actually S(τ ) is independent of τ and is equal to the instanton action S inst , Eq. (B4).Therefore the GR-QTST rate will not depend on the choice of τ , and will be given by which is also the instanton and exact Fermi's goldenrule result for a system of two linear potentials. 46In the classical limit, the prefactor is the same, but the exponent disappears as S inst → 0. Interestingly if we had ignored the fluctuation terms, GR-QTST would not have given the correct result.This shows that it is not just the instanton, but also the fluctuations around all the constrained minima which contribute to the path integral and to the success of the new method.
This proof is valid not just for systems with κ 0 κ 1 < 0, but also for those with κ 0 κ 1 > 0, which describe the Marcus inverted regime.Note that the ring-polymer instanton method is not stable numerically in the inverted regime, whereas GR-QTST is.
Although the reactant partition function, Z 0 , is not well defined for this system, this important model system describes a limiting case of many nonadiabatic reactions, for which the reactant partition function can be defined.Our findings here that GR-QTST gives excellent predictions for the rate constant are expected to carry over at least approximately to more general systems.
Appendix C: Analytic solution for an uncoupled linear and harmonic system In this section, we consider a multidimensional system formed of crossed linear potentials uncoupled to harmonic oscillators.We choose it to be uncoupled such that the analytical result can be obtained.The numerical result of a coupled problem is given in Sec.V B.
Here we just consider the effective action of the y-part, which is defined as follows: e −φ (y) (τ )/ = J 0 Λ −N0 0 where where c = mω 2 , which is just a constant independent of index k.
The constraint functional for the uncoupled linear and harmonic potential-energy surface is the sum of the energy-matching functionals in x and y directions, σ (x,y) (τ ) = σ (x) (τ ) + σ (y) (τ ). (C7) where σ (x) (τ ) is the constraint functional for the onedimensional problem, Eq. (B14).For the end points of the classical path, the energy matching function for the harmonic oscillator simply is σ (y) (τ ) = 2 3 β η (C12) The integral over end points, y and y , are also Gaussian, and yields the following result: (C15) where θ n = csch(ωτ n ) (ω coth(ωτ n ) − 1/τ n ).Therefore the effective action is defined by e −φ (y) (τ )/ = m τ 0 τ 1 Under the GR-QTST formulation, the effective action is evaluated at τ = τ * , which corresponds to the point where φ (y) (τ ) demonstrates the maximum.For this system, this maximum occurs at τ * = τ inst for which χ τinst = 0. Therefore the effective action simplifies to e −φy(τ * )/ = 1 2 sinh(β ω/2) . (C17) This quantity is recognized as the partition function in the y-direction, which cancels with an equivalent term in the reactant partition function.Therefore the GR-QTST rate obtained for an uncoupled linear and harmonic system is exactly the same as the rate obtained for a one-dimensional linear system, which is of course correct.This derivation also holds for many uncoupled oscillators and remains at least approximately true for a coupled bath (see Fig. 5).

8 FIG. 1 .
FIG.1.Density plots of the relative free energy, Fτ (E) = − log pτ (E), over a range of τ values.pτ (E) is the probability density of finding the average hopping point, x+, with the polarization E ≡ 1 2 V−(x+) at a given τ obtained from path-integral sampling with (a) no constraint, (b) the thermodynamic constraint, and (c) the virial constraint.(Note that here the probability is normalized for each value of τ such that pτ (E) dE = 1.)This system is the one-dimensional harmonic model defined in Sec.V B, with /λ = 0.5 and no bath modes.For this system, τinst/β ≈ 0.79.The regions which are rarely sampled are coloured yellow whereas the commonly sampled regions are blue.The crosses signify the instanton configurations.The insets show the path with the minimal action at the indicated values of τ and E. The x-axis of the insets is the nuclear configuration, x, and the y-axis shows the potential-energy surfaces and the instantaneous energy of the beads.

)
FIG. 5.The dependence of the effective action on τ for the system-bath model with /λ = 0.5 and in each case, a different number of dimensions, f .All calculations were performed using N = 50.GR-QTST results, computed using Eq.(32), are shown with statistical error bars.The result of the steepest-descent (SD) approximation, Eq. (35), is included only for f = 8 as for all multidimensional (f ≥ 2) cases, the SD approximations to φ(τ ) are approximately the same.The effective actions of the instanton method, φinst, are marked with dots at τ = τinst for each case, although these lie almost exactly on top of each other.All results are normalized by the partition function Z0 to ease comparison with the frictionless results shown in Fig.3.

FIG. 6 .
FIG.6.The calculated rate constants are shown as a function of bias, , for the 8-dimensional system-bath model.The various approaches used are exact Fermi's golden rule,43 classical Marcus theory, Eq. (51), semiclassical instanton (SC), Eq. (52), and GR-QTST.Wolynes theory is equal to SC for this case but is only numerically applicable in the normal regime, < λ.The reference, kMT(0), is the classical Marcus theory rate constant for = 0. SC rates were computed analytically using N → ∞, whereas the GR-QTST rates were obtained with N = 100 ring-polymer beads.