Semiclassical instanton formulation of Marcus-Levich-Jortner theory

Marcus-Levich-Jortner (MLJ) theory is one of the most commonly used methods for including nuclear quantum effects into the calculation of electron-transfer rates and for interpreting experimental data. It divides the molecular problem into a subsystem treated quantum-mechanically by Fermi's golden rule and a solvent bath treated by classical Marcus theory. As an extension of this idea, we here present a"reduced"semiclassical instanton theory, which is a multiscale method for simulating quantum tunnelling of the subsystem in molecular detail in the presence of a harmonic bath. We demonstrate that instanton theory is typically significantly more accurate than the cumulant expansion or the semiclassical Franck-Condon sum, which can give orders-of-magnitude errors and in general do not obey detailed balance. As opposed to MLJ theory, which is based on wavefunctions, instanton theory is based on path integrals and thus does not require solutions of the Schr\"odinger equation, nor even global knowledge of the ground- and excited-state potentials within the subsystem. It can thus be efficiently applied to complex, anharmonic multidimensional subsystems without making further approximations. In addition to predicting accurate rates, instanton theory gives a high level of insight into the reaction mechanism by locating the dominant tunnelling pathway as well as providing information on the reactant and product vibrational states involved in the reaction and the activation energy in the bath similarly to what would be found with MLJ theory.


I. INTRODUCTION
The theoretical study of electron-transfer is of essential importance and relevance not only because these reactions are a key step in many chemical and biological processes but also because the methods developed to deal with them can be applied in many other scenarios ranging far beyond their original scope. This follows from the fact that electron-transfer reactions are just one example of the more general set of curve-crossing problems. Hence, contributions to the understanding of electrontransfer reactions have been made with various motivations including electrochemistry, molecular spectroscopy, polaron transport as well as more general atom-transfer reactions, which led to different ways of tackling the problem from classical dielectric continuum theory to a full quantum molecular picture. 1, 2 Inspired by earlier work, 3 Marcus based his theory of electron transfer, for which he later won the Nobel prize in 1992, 4 first in terms of a dielectric solvent continuum 5 and later on a classical statistical mechanical description of the solvent. 6 To this day, Marcus theory is probably the most commonly applied approach for the description of electron-transfer reactions and initiated tremendous development involving electron and hole transfer between atoms, molecules or even proteins, in the condensed phase as well as at interfaces. 7,8 Hence, his findings had and still have an enormous impact on a multitude of scientific disciplines comprising solution chemistry, solid-state physics as well as biological processes. 9 One of the essential insights from Marcus' classical theory was the prediction of the so called "inverted a) Electronic mail: eric.heller@phys.chem.ethz.ch b) Electronic mail: jeremy.richardson@phys.chem.ethz.ch regime", 6 the existence of which was later confirmed by experiment, 10 where the rate decreases as the thermodynamic driving force grows larger than the reorganization energy. Soon, however, it was realized by theory and experiment that the neglect of nuclear quantum effects in Marcus theory can lead to dramatic errors of several orders of magnitude in the rate, especially in the inverted regime. [11][12][13] Based on the connection to spectroscopy and solidstate nonradiative processes Levich and coworkers put the theory onto a rigorous quantum-mechanical basis and introduced a quantum statistical mechanical description of outer sphere electron transfer. 14 This was done by employing Fermi's golden rule 15 formula for the quantum transition rate, which is obtained as the nonadiabatic (weak-coupling) limit from perturbation theory. 16,17 However, because outer-sphere electrontransfer is typically dominated by the low-frequency solvent modes, the resulting quantum effects are rather small.
Several years later, crucial advancements were made in particular by Jortner and coworkers by explicitly taking the reorganization of the inner sphere into account. 1,2, [18][19][20][21][22] As opposed to the solvent, the inner sphere often exhibits intra-and intermolecular rearrangements associated with high-frequency vibrational modes which are therefore subject to substantial quantum effects. Hence, they treated the inner sphere quantummechanically using Fermi's golden rule while keeping the classical approximation for the solvent bath. The resulting Marcus-Levich-Jortner (MLJ) theory constituted a considerable progress for the whole field, as it was the first rigorously derived method able to describe nuclear quantum effects in electron-transfer reactions which was valid throughout virtually the whole temperature range. 23 Thus, the method poses a vital step towards the goal of establishing a unified description of electron transfer in the various fields mentioned above across diverse time and temperature scales. 8,24,25 MLJ theory is broadly applied to the prediction and explanation of charge-carrier mobilities 26 often with the objective to give a guideline for the synthetic study and reasonable design of high-performance semiconductors that can be applied in organic photovoltaics. [27][28][29][30] The application to large systems can be facilitated by using the theory in conjunction with density-functional theory. 31 Furthermore it can be applied in the study of molecular junctions, 32 photonics, 33 polaritons 34 and polarons 35 as well as for the description of spin transitions and phosphorescence. [36][37][38][39][40] Besides these technological disciplines, it is also frequently applied for the understanding of complex chemistry, 41 electron transfer in supermolecules 42 and biochemistry, [43][44][45] charge transfer in DNA [46][47][48] and photosynthesis. 49,50 As tunnelling is especially prevalent in the Marcus inverted regime, 51 MLJ theory is of particular interest in the study of molecular electron-transfer reactions which are strongly exothermic [52][53][54][55] or which are initiated by photoexcitation. 56,57 The MLJ description of quantum tunnelling, which will be extensively discussed in the later sections, is understood by shifting the Marcus parabolas (free energy along the bath coordinates) by the quantized energy levels of the inner sphere. The rate therefore consists of contributions from multiple vibrational channels weighted according to a thermal distribution. 58 This interpretation appears quite different from the standard picture of tunnelling in which a particle penetrates a potential energy barrier with an energy smaller than the barrier height. The main disadvantage of the approach is that it requires wavefunction solutions of the time-independent Schrödinger equation in order to compute the energy levels and Franck-Condon overlaps, which severely limits its usefulness for the description of realistic, anharmonic and multidimensional systems.
A number of alternative methods which can be used for the study of electron-transfer reactions and other goldenrule processes [59][60][61][62][63][64][65][66][67][68][69][70][71][72][73][74][75] are based on Feynman's path-integral description of quantum mechanics 76,77 rather than on wave mechanics. From this set, semiclassical golden-rule instanton theory 78,79 in particular bears multiple appealing features. Without any prior knowledge about the analytic shape of the potential, it locates the "instanton" in the full-dimensional configuration space of the system, which can be thought of as the optimal tunneling pathway, [80][81][82][83] and therefore provides direct insight into the reaction mechanism. Furthermore the method was recently extended towards the Marcus inverted regime, 84 which otherwise typically poses a problem for imaginary-time path-integral approaches, 85 although some extrapolation techniques have been used successfully to avoid this problem in other methods. 67 By employing a ring-polymer discretization to the paths, 86 the instanton method is able to simulate tunnelling in multidimensional, anharmonic systems in a computationally efficient way and is ideally suited for calculations in conjunction with high-level electronic structure methods just as in the standard adiabatic formulation of the theory. 81,[87][88][89][90][91][92][93][94] Golden-rule instanton theory constitutes a semiclassical path-integral formulation of Fermi's golden rule and hence has the potential to be applied in a multitude of different fields, just as the golden rule itself.
One of the great strengths of instanton theory is its full-dimensional formulation of tunnelling such that it does not rely on an a priori choice of the reaction coordinate. However, for many relevant reactions, especially in the condensed phase, even if one does not know the exact tunnelling path, one already has a good idea of which part of the system under investigation has to be considered explicitly and which part can be accounted for on a coarser level. It is this same separation into an inner and outer sphere which was the cornerstone of MLJ theory. Therefore in this paper, the formalism for a "reduced" semiclassical golden-rule instanton theory will be laid out, which describes tunnelling within the modes of the inner sphere under the implicit influence of either a classical or quantum harmonic bath. The presence of the bath affects the equations of motion of the inner sphere and renders the resulting reduced instanton non-energyconserving due to energy exchange between inner and outer sphere. This is analogous to the reduced density matrix formalism employed in the study of open quantum systems. The resulting instanton picture preserves the convenient interpretation of quantum tunnelling as a particle travelling in the classically forbidden region below the barrier.
Although the formalisms seem at first glance rather different, we will draw a connection between the MLJ and instanton theories by deriving them both from a common expression. In doing so, it will be shown clearly that the instanton approximation is fundamentally different from other approximations such as the broadly applied cumulant expansion method 95 and the semiclassical Franck-Condon sum. 11,96 Numerical results demonstrate that instanton theory is very accurate over a range of systems including anharmonic modes where these alternative approximations break down.
Some rate theories have the advantage that they are based on expressions which are simple enough that one can easily see their dependence on certain parameters and thus gain insight into the behaviour of different systems. 2, 97 We will argue that instanton theory allows for a well-balanced combination of easily attainable insights, as well providing a realistic molecular simulation. Even when applied to complex anharmonic multidimensional potentials, the method uniquely identifies an optimal tunnelling pathway which provides a simple onedimensional picture of the reaction, highlighting which modes are involved in the tunnelling. In addition to this, the instanton can be analysed to obtain information on the energies of the initial and final states of the system before and after the electron-transfer event, similar to what is computed in MLJ theory as we will show.

II. GOLDEN-RULE RATE
The total Hamiltonian which describes electron transfer between a reactant |0 and product |1 electronic state is defined by 98 where the electronic interaction between the states is given by the nonadiabatic coupling ∆. Throughout this work the electronic coupling is taken to be constant, but the generalization to position-dependent couplings is fairly straight forward. 78,99 Furthermore the coupling is assumed to be very weak such that the rates occur in the golden-rule limit, i.e. ∆ → 0, which is typically the case in electron-transfer reactions. 98 A driving force, ε, has been included explicitly in the total Hamiltonian, which could describe an internal energy bias or the effect of an external field. It is kept separate here for clarity but it could of course be simply absorbed into the definition of H 1 . We are interested in studying problems which can be subdivided into an inner sphere, whose molecular structural characteristics will be explicitly taken into account, and an outer sphere, which typically includes the solvent degrees of freedom and will be treated as an effective harmonic environment characterized by its spectral density. In the language of open quantum systems, these are called subsystem and bath and are here taken to be uncoupled to each other, 100 although there is of course still some coupling through the nonadiabatic terms in Eq. (1). Hence, the full nuclear Hamiltonian for electronic state |n can be written asĤ The subsystem Hamiltonians,Ĥ s n , only depend on the coordinates q = (q 1 , . . . , q d ) and their conjugate momenta p = (p 1 , . . . , p d ), while the bath Hamiltonians,Ĥ b n , are solely a function of the coordinates Q = (Q 1 , . . . , Q D ) and momenta P = (P 1 , . . . , P D ). Without loss of generality, these degrees of freedom have been mass-weighted such that all subsystem modes are associated with the same mass, m, and likewise all bath modes with mass M .
The harmonic approximation for the bath will be employed: 1 where the plus sign corresponds to the reactant state and minus sign to the product state. The bath Hamiltonians thus combine in Eq.
(1) to describe a spin-boson model, 101 defined by the associated frequencies {Ω j } and displacements {ζ j }, which can be selected such that they represent an appropriate spectral density. This spinboson model is complemented by the subsystem modes, whose potential-energy surfaces will be kept general for the derivations in this work such that they can, in principle, provide a realistic description of an anharmonic molecule.
The full system is prepared as a thermal equilibrium ensemble in the reactant state with inverse temperature β = 1/k B T and partition function Z 0 = Tr e −βĤ0 . The quantum-mechanical rate expression for a reaction from the reactant to the product electronic state in the goldenrule regime can be derived from a perturbation expansion to lowest order in the nonadiabatic coupling ∆ between the two electronic states of an integral over the flux correlation function 98,102 to give The flux-correlation function is an analytic function of time, and hence, the rate is independent of the imaginary-time parameter τ , 103 although it has been included explicitly as it will play a pivotal role in the semiclassical approximations taken later on. Due to separability of the subsystem and bath, a quantum trace can be taken independently over the respective contributions. The reactant partition function thus factorizes according to Z 0 = Z s 0 Z b 0 into a subsystem part, Z s 0 = Tr s e −βĤ s 0 , and a bath part, Z b 0 = Tr b e −βĤ b 0 . The correlation function likewise splits into product of subsystem and bath parts.
The rate can thus be rewritten using the convolution theorem of Fourier transforms: 20 where the subsystem and bath lineshape functions are The equivalence to Eq. (5) can easily be checked by substituting Eqs. (7) into Eq. (6) after renaming the integration variable t in the two cases to t s or t b and using The lineshape function of the subsystem expanded simultaneously in the position and eigenstate bases can be written as where E µ 0 and E ν 1 are the internal energy levels and ψ µ 0 and ψ ν 1 are the corresponding wavefunctions of reactants and products, respectively.
Because of the global harmonic approximation of the bath, the well-known result for the spin-boson model 1, 101 can be used to cast Eq. (7b) into where the effective action of the bath is defined by Note that this is an analytic function of its argument and can therefore also be used to describe real-time dynamics in Eq. (9). The coordinate dependence of the bath has been completely integrated out, which is the reason why the effective action [Eq. (10)] only depends on time. Assuming the spectral density of the bath is known, the time-integral of Eq. (9) can be carried out (either by quadrature or by steepest descent) in order to account for quantum effects within the solvent. 60,104 In cases where the bath represents a polar solvent environment, which typically comprises long-wavelength polarization modes, it is often justified to approximate the effective action [Eq. (10)] by its classical, low-frequency limit where |Ω j τ | 1 and β Ω j 1. 1 In this case, the classical bath action is where the bath reorganization energy is given by In these formulas, τ /β plays the role of a "symmetry factor" as described in Ref.

2.
In order to include a quantum harmonic bath with Eq. (10), knowledge of the bath spectral density is required to define {Ω j } and {ζ j }. On the other hand, a classical harmonic bath [Eq. (11)] can be simpler to employ as it is fully characterized by its reorganization energy Λ b and thus requires much less information.
The approach which we will follow in this paper is to evaluate the subsystem and bath lineshape functions using different representations and approximations to derive multiple methods for computing electron-transfer rates in the golden-rule regime.
For instance, if the relaxation of the inner sphere is assumed to play no role in the reaction under consideration, there is no subsystem contribution to the Hamiltonians in Eq. (2). The subsystem lineshape function Eq. (7a) therefore reduces to I s (v) = e +(τ +it)v/ dt = 2π δ(v). Employing the classical approximation for the action [Eq. (11)] in the bath lineshape function Eq. (9), plugging the lineshape functions into Eq. (6) and performing the final time-integral analytically leads the famous Marcus rate equation 13 It describes electron-transfer reactions which do not involve significant rearrangements within the inner sphere (subsystem) and are therefore determined only by the conformational changes in the outer sphere (bath). As is well known, Marcus theory thus gives the correct classical limit of the rate in the case of a spin-boson model. 1, 14,105 In an alternative and more powerful derivation of Marcus theory, the trace could have been evaluated over the bath degrees of freedom in Eq. (7b) directly by a classical phase-space integral. 85,106 In fact we could treat the subsystem in the same way to obtain k MT (ε), a theory equivalent to Eq. (13) but written in terms of the total reorganization energy, Λ = Λ s + Λ b . This treatment allows for anharmonic potential-energy surfaces but reduces to give the same rate formula as long as the free-energy surfaces themselves are harmonic. In this case, one should treat the driving force ε as a free energy as it can also include entropic effects. 107 In many cases, however, the inner sphere undergoes significant conformational changes as well and can therefore not be ignored. Moreover, the molecules in the reaction center commonly exhibit high-frequency modes, which necessitates the explicit consideration of quantum effects (such as the existence of zero-point energy and possibility of tunnelling) within an anharmonic environment. We can derive various methods to compute the subsystem contribution simply by carrying out the sums and integrals of Eq. (8) in different orders. Although all these approaches give identical results in their exact form, they provide different starting points for taking approximations.

III. MARCUS-LEVICH-JORTNER THEORY
The subdivision of the full nuclear Hamiltonians of each electronic state into independent subsystem and bath parts [Eq. (2)] is the foundation on which MLJ theory is grounded. 18,19 In this approach one then treats the subsystem quantum mechanically and the bath classically.

A. Formalism
Here we rederive MLJ theory on the basis of the formalism laid out in Sec II. Starting from Eq. (8), we first take the integrals over positions and time and set τ to be zero. Consequentially one arrives at Fermi's goldenrule (FGR) formula 108 for the lineshape function of the subsystem 19,20 where θ µν = ψ µ 0 (q) * ψ ν 1 (q) dq are the Franck-Condon factors and the subsystem contribution to the reactant partition function in the energy eigenbasis is Z s 0 = µ e −βE µ 0 . Combining this wavefunction representation of the subsystem part with the bath lineshape function [Eq. (9)] using the classical effective action [Eq. (11)] and performing the final convolution integral in Eq. (6) leads directly to the Marcus-Levich-Jortner electron-transfer rate theory in a system of two crossing potentials of arbitrary shape in a classical harmonic bath 20 with which is the most general version of MLJ theory used in this work. The total rate in Eq. (15) comprises contributions from all reactant and product vibrational channels. This "static" formulation of electron transfer (i.e. time has been integrated out) results in a rate expression that requires knowledge of all internal states of the subsystem Hamiltonians,Ĥ s n . For complex anharmonic molecules, this is not possible to compute without further approximations. Thus, the most commonly employed form of the Marcus-Levich-Jortner theory takes the extra approximation that the subsystem potentials for the reactant and product are displaced one-dimensional harmonic oscillators with identical frequencies, ω. Motivated by the fact that in many problems of physical interest the subsystem comprises very high-frequency modes, it is often appropriate to assume that the thermal energy is low compared with the energy spacing in these modes. Hence, only transitions from the ground vibrational reactant state with quantum number µ = 0 have to be considered. The general expression for the onedimensional overlap integral of two displaced harmonic oscillator wavefunctions can therefore be further simplified, because only the terms 22 with A = Λ s / ω, have to be taken into account. This results in the well known rate formula for a single quantum harmonic mode in the low-temperature limit 21 where the rate into the product-state ν is This low-temperature rate therefore consists of contributions from multiple, parallel product vibrational channels with effective driving forces of ε ν = ε − ν ω. As can be seen from the exponential "activation" part of Eq. (19), the major contributions to the rate will typically involve the product vibrational states whose effective driving force ε ν is approximately equal to the bath reorganization energy. In cases where ε < Λ b , this is not possible, and so then the ν = 0 product state is expected to dominate. Thus, the dominant product vibrational state will depend on the thermodynamic driving force and transitions to highly excited vibrational states are of particular importance for very exothermic reactions and hence especially in the inverted regime.

B. Model example
The simple model which we will use to illustrate MLJ theory is formed of two-dimensional displaced harmonic oscillators with one mode treated quantum mechanically and the other classically. The subsystem potentials are defined by where the reactant state is associated with the plus sign and the product state with the minus sign, and because In particular we will apply the theory to two different models, defined by the parameters in Table I, one of which is in the normal and the other in the inverted regime.
The parameters are chosen so as to illustrate two common scenarios. Because the inner sphere typically TABLE I: Definition of the parameters used in the two harmonic models studied in this work. The rate is independent of the masses m and M , which therefore do not have to be defined. As is common practice, the value of the frequencies are defined by their related wavenumber.
1000 500 Ω (cm −1 ) 50 50 comprises high-frequency vibrational modes, the lowtemperature limit of the MLJ rate [Eq. (18)] can often be applied to a good approximation. 1,20 Hence, activated vibrational reaction channels only have to be considered for the product. This case is exemplified by the invertedregime model. The situation is illustrated in Fig. 1(a) which shows the product potential shifted by the vibrational energy gap ν ω. In Fig. 1(c) the reaction rate is broken down into contributions from the individual product channels, which are clearly centered around the dominant vibrational state ν = 15 and rapidly fall off on either side. Sometimes, however, a system also requires the consideration of activated vibrational states of the reactant, which necessitates the use of the general expression Eq. (15). As illustrated in Fig. 1(d), this is the case for the normal-regime model, where excited vibrational states of both the reactant and product make significant contributions to the rate. Thus, in Fig. 1(b), not only the product but also the reactant potential is shifted by the vibrational energies. The required Franck-Condon overlap integrals for a subsystem of two displaced harmonic oscillators can be computed with well-known analytic formulas. 95 Fig. 1(d) shows that the dominant contribution to the rate comes from the reaction channel from µ = 3 to ν = 11.
Perhaps even more important than the ability of MLJ theory to predict rates is that it provides this simple picture of the quantum nuclear effect on an electrontransfer reaction. By viewing the reaction along the bath coordinates and shifting the potential-energy surfaces by the excitation energies of the subsystem, one obtains vibrational-state resolved contributions to the rate, which are centred around a dominant vibrational channel. This is a popular way of understanding reactions and has for instance been used to explain why rates in the inverted regime commonly flatten off instead of decreasing rapidly with driving force as predicted by classical Marcus theory [Eq. (13)]. 10,12 In the inverted regime, it can easily be seen from Eq. (18) and Figs. 1(a) and (c) that the dominant contribution to the rate originates from the vibrational channel that approximately shifts the bath product potential to the activationless regime, where ε ν ≈ Λ b , 109 and thus predicts a rate approximately independent of driving force. 110 This analysis of the inverted-regime model is based on the simplification that the rate is fully determined by the exponential "activation" part. In reality, however, rates in the inverted regime are also affected by the Franck-Condon factors such that they do not actually become constant with driving force. We find that the bath activation energy for the dominant vibrational transition in the inverted-regime model is 0.51 kcal mol −1 , which is almost activationless, but still not negligible relative to the thermal energy. In the normal-regime model, where excited reactant states also play an essential role as illustrated in Figs. 1(b) and (d), the full rate expression Eq. (15) is no longer dominated by an activationless channel at all. We find a significant activation energy for the dominant vibrational channel of 6.64 kcal mol −1 , which illustrates the compromise between minimization of the activation energy and maximization of the Franck-Condon overlaps that has to be made. This considerably complicates the interpretation of the MLJ rate formula even when the harmonic oscillator approximation is employed.
For more realistic systems described by multidimensional anharmonic potential-energy surfaces the Franck-Condon factors are practically impossible to obtain and the subsystem is thus commonly approximated by simple models for which these are known analytically. This introduces unknown errors into the predicted rate, and it is to avoid this problem that we now turn to instanton theory.

IV. REDUCED INSTANTON THEORY
Inspired by the Marcus-Levich-Jortner approach we will derive a reduced instanton theory, where only the inner sphere is treated explicitly in molecular detail while the outer solvent shells are accounted for with the harmonic bath approximation.

A. Formalism
In order to derive the semiclassical golden-rule instanton rate expression, we start from the time-dependent correlation function formulation of the reaction rate in Eq. (5). The trace can be split up into a subsystem and bath contribution, where the latter can, due to its harmonic nature, again be replaced with the well known solution for the spin-boson model in terms of the effective bath action Φ(τ ). In contrast to the MLJ approach, the trace in the subsystem coordinates will be expanded in the position basis, which leads the following expression for the rate: Again in analogy to the dynamics of open quantum systems, e −Φ(τ +it)/ plays the role of an influence function. 77,101 The matrix elements of the quantum propagators, describe the dynamics of the subsystem variables evolving according to the HamiltoniansĤ s n from the initial positions q i to the respective final position q f in imaginary time τ n . The imaginary-time propagators are equivalent to quantum Boltzmann distributions and it is this connection which allows instanton theory to approximate the thermal rate in a statistical way using imaginary-time dynamics.
If the imaginary-time propagators and spatial integrals were evaluated by path-integral Monte Carlo calculations and the remaining time integral taken by steepest descent, one would obtain a version of Wolynes theory where the bath is treated implicitly by the influence function. 59,68 This, however, is not the purpose of this work as we wish to derive a semiclassical instanton formulation of the rate.
Instead we replace the quantum propagators by the corresponding van-Vleck propagators 111 generalized for imaginary-time arguments 87,112 thus introducing a semiclassical approximation. The resulting expression is evaluated by locating the classical trajectory, q n (u), travelling in imaginary time u, which makes the Euclidean action of the subsystem, S n , stationary. The action for a path travelling from its initial position q n (0) = q i to its final position q n (τ n ) = q f in imaginary time τ n is defined as whereq n (u) = dqn du is the imaginary-time velocity. The prefactor of the semiclassical propagator is given by the determinant By multiplying the two propagators in Eq. (21) together, we obtain the total action S(q , q , τ ) = S 0 (q , q , β − τ ) + S 1 (q , q , τ ), (26) as the sum of contributions from two trajectories, one of which travels on the reactant potential and the other on the product potential. These trajectories join each other to form a continuous periodic pathway, called the instanton. The imaginary times τ n associated with the two paths are given by τ 0 = β − τ and τ 1 = τ .
Combining this result with the effective action of the bath according to Eq. (21), the total effective action becomes where one could employ either the effective quantum bath action from Eq. (10) or its classical limit Eq. (11). The only effect of the bath is to thus alter the total action by adding an extra τ -dependence alongside the driving force term. However, as we will show, the simple addition of the bath action can lead to significant changes for the instanton path and for our interpretation of the reaction mechanism.
In order to obtain the semiclassical instanton expression for the rate, the integrals over q and q as well as the time-integral will be carried out by steepest descent. Therefore it is necessary first to study the path corresponding to the stationary point of the effective action, for which ∂ S r ∂q = ∂ S r ∂q = ∂ S r ∂τ = 0. This path is our definition of the reduced instanton and by analyzing the consequences of vanishing derivatives, we can understand its properties in the general case.
As in the standard golden-rule instanton formulation, 78 the instanton pathway consists of two trajectories q 0 (u) and q 1 (u), which join smoothly into each other at the stationary or hopping point in the subsystem coordinate space q = q = q ‡ . Because the q-derivatives are not altered by the influence of the bath, the momentum, given by p = − ∂ S0 ∂q = ∂ S1 ∂q (equivalent for double primes), is thus still conserved across the hopping point. In this sense, the instanton therefore remains a periodic orbit of imaginary time β as in the standard theory. 78 However, a major difference occurs due to the bath's influence on the derivative with respect to τ . The subsystem energies of the two trajectories are given by Therefore, in the case without the presence of a bath, the condition at the stationary point is given by Considering ε as a contribution to the product energy as was done in Ref. 78, this relationship implies that the reaction conserves energy, i.e. E s 0 = E s 1 − ε. The hopping point must therefore be located on the crossing seam where V 0 (q) = V 1 (q) − ε.
This no longer holds true once bath modes are added. Then the condition at the stationary point changes to The presence of the bath will thus affect the stationary value of τ and hence the entire instanton path and the value of its action. 113 In particular, the energies of the two trajectories no longer match E s 0 = E s 1 − ε in general. Hence, the presence of the bath renders the reduced instanton non energy-conserving within the subsystem. However, the energy change in the subsystem is exactly compensated by an energy change of opposite sign in the bath given by ∆E b = ∂ Φ ∂τ such that ∆E s + ∆E b − ε = 0. This is in agreement with what one would expect from an open quantum system, in which only the total combined energy of subsystem and bath is conserved but not the individual components. In this theory one does not have direct access to the energies of the bath which would be needed to fully justify this interpretation for ∂ Φ ∂τ . However, we will show that this definition is correct in Sec. V B.
Potential curves V s 0 (q) (blue solid lines) and V s 1 (q) − ε (orange solid lines) of the subsystem as defined in Eq. (32) with ε = Λ/4 = 31.7 kcal mol −1 . The two trajectories of the reduced instanton are shown at discrete time steps by blue (reactant) and red (product) dots, and their energies, E s 0 and E s 1 − ε, are depicted by the blue and red dashed lines. These energies are separated by the energy gap ∆E s − ε, which is equal to the potential-energy gap at the hopping point (q ‡ , purple dot). The inset shows the instanton for a model with the same subsystem but without the presence of a bath enlarged from the area framed by the grey dotted box. Here ∆E s − ε = 0 and therefore energy conservation is satisfied and the hopping point (purple dot) is located where the potentials cross.
One consequence of the energy jump caused by the presence of the bath is that the hopping point, q ‡ , is not located on the crossing seam between the two subsystem potentials. In fact, because the momenta of the two trajectories are equal at the hopping point, the energy jump within the subsystem must correspond exactly to the potential energy difference, ∆E s = ∆V s (q ‡ ), where ∆V s (q) = V s 1 (q) − V s 0 (q). In Fig. 2, we illustrate the reduced instanton pathway for the anharmonic model discussed in Sec. IV B, as well as the the energies of the two trajectories as defined by Eq. (28).
Note that the reactant or product energy is conserved along its respective trajectory and is thus identical to the potential at the turning point, which can be easily seen in the figure as the point with lowest potential along the path. The concept of a turning point in this context can be understood by the fact that dynamics in imaginary time are equivalent to real-time dynamics on the upside-down potential. 112 At the turning point, the paths therefore bounce against the potential which they are travelling on.
Because the instanton orbit folds back on itself and is therefore not so easy to depict, it is worth describing it in a little more detail. If we first follow the instanton path-way in Fig. 2 along the trajectory q 0 (u) starting at the (purple) hopping point on V s 0 and with a certain amount of momentum pointing to the left, we find that the path descends towards the reactant state minimum, where it bounces against the potential and returns to where it started but with momentum now pointing to the right. So far this is equivalent to the instanton pathway in the standard formulation of the theory. 78 Once the hopping point is reached, however, a sudden jump in potential energy occurs, which accompanies the transition into the product state. This is in stark contrast to the standard formulation, where both trajectories tunnel at the same energy, as shown in the inset of Fig. 2 for a case without a bath. After the state transition, we travel along the trajectory q 1 (u), which after reaching the turning point on V s 1 also returns to the hopping point. A transition back to the initial hopping point on V s 0 completes the periodic cycle.
After the instanton pathway has been located, the integrals in Eq. (21), where the propagators have been replaced with Eq. (23), can be carried out by steepestdescent integration around the stationary point. Thus we arrive at the reduced instanton expression for the goldenrule rate where all quantities are evaluated at the stationary point of S r (q , q , τ ) except the reactant partition function Z s 0 , which is treated by an equivalent steepest-descent approximation around the minimum of the reactant. 87 The additional prefactor from the steepest descent integration in the subsystem positions evaluates to the determinant where we have used the fact that derivatives of S r with respect to the end points are equal to derivatives of S.
The bath therefore has no direct effect on C but does explicitly appear in d 2 S r dτ 2 = d 2 S dτ 2 + d 2 Φ dτ 2 as well having an important effect on the instanton path itself as previously discussed. Apart from these changes, the formula resembles the semiclassical golden-rule instanton rate expression derived in previous work 78 and gives identical results without needing to treat the harmonic bath explicitly.
Just as for previous golden-rule instanton calculations, 79,84 a ring-polymer discretization scheme of the instanton pathway is employed in order to describe nonadiabatic reactions for multidimensional, anharmonic systems. By adopting the ring-polymer formalism, the localization of the instanton path, which is defined as a stationary point of the action in Eq. (27) in the coordinate and τ variables together, reduces to a standard saddle-point search problem which can be solved numerically with well-established optimization algorithms. Algorithms for computing the necessary derivatives of the action as well as detailed information about the optimization scheme can be found in Ref. 86.
In our recent extension of the theory, 84 we have shown that ring-polymer instanton theory can equivalently be utilized to compute electron-transfer rates in the Marcus inverted regime, where tunnelling effects commonly play a particularly important role. The major difference in this regime is, that one of the two paths travels in negative imaginary time, which allows an analogy to the physics of antiparticles. 114 In the computational realization, this difference manifests itself merely in a slight change of the optimization algorithm. Hence, whereas in the normal regime the instanton is a single-index saddle point of the ring-polymer action in the combined space of ring-polymer coordinates and imaginary time, in the inverted regime the instanton path corresponds to a higherindex saddle point of the ring-polymer action. The index of a saddle point here defines the number of negative eigenvalues in the second-derivative matrix of the ringpolymer hessian at this point. But since we exactly know the index of the desired saddle-point, the instanton can be optimized with the same routines by using standard eigenvector-following schemes. We thus take uphill steps in the direction of eigenvectors corresponding to negative eigenvalues and standard down-hill steps in the direction of eigenvectors associated with positive eigenvalues. This methodology can be directly transferred to the reduced instanton picture without any additional complications and hence allows us to apply it to the normal and inverted regimes alike.
The advantage of the reduced instanton approach is that the optimization is confined to the inner sphere and τ -coordinates only, whereas the only direct influence of the bath on the optimization procedure manifests itself in an external field in the imaginary-time variable. This reduces the computational costs of the simulation and enables it to be applied within a multiscale modelling approach, where certain parts of a system are treated at higher levels of accuracy than others.

B. Model example
We will employ the newly formulated reduced instanton method along with MLJ theory to compute reaction rates of an anharmonic subsystem of two bound Morse oscillators in a multidimensional harmonic bath. The subsystem is defined by the potentials (depicted in Fig. 2) where α 0 = 1.5Å −1 and α 1 = 1.4Å −1 determine the length scales, ξ 0 = 1.0Å and ξ 1 = 1.5Å are the equilibrium positions, and D e 0 = 115 kcal mol −1 and D e 1 = 80 kcal mol −1 are the dissociation energies of reactants and products. The (product) reorganization energy of this subsystem is therefore Λ s = V 1 (q min is the minimum of V n (q). The reduced mass is chosen to be m = 1.10 u. The well frequencies of the two Morse oscillators obtained by harmonic analysis are ω n = α n 2D e n /m, which results in frequencies of ω 0 = 2358 cm −1 and ω 1 = 1835 cm −1 for the reactant and product well respectively. The Schrödinger equation for the Morse oscillator can be solved analytically to give the bound-state energies and likewise for E ν 1 , where the dimensionless anharmonicity parameters of the Morse oscillators are defined by χ n = α 2 n /2mω n and in this case have values of 0.015 and 0.016 for the reactant and product potential respectively.
The bath is defined by the discretized spectral density and the D = 100 bath modes were chosen according to Ref. 115 as The effective mass of the bath modes M does not have to be specified, as the rate is independent of this choice. The frequency spectrum is bounded from above by the maximum frequency Ω max = 3000 cm −1 and thus has the density 116 The couplings were chosen to emulate a Debye spectral density defined by with characteristic frequency Ω c = 500 cm −1 and η = 25 kcal mol −1 . Hence the coupling constants c j are determined by the formula which are related to the shifts ζ j in Eq. (4) by c j = M Ω 2 j ζ j . The reorganization energy of the bath is then obtained by Eq. (12), which in our case results in Λ b = 44.8 kcal mol −1 . The total reorganization energy of the subsystem and bath combined is therefore given by Λ = Λ s +Λ b = 127.0 kcal mol −1 . The temperature for the rate calculations was chosen to be 300 K.
Similar models were studied with MLJ theory in Ref. 117. In order to make use of analytical formulas for the Franck-Condon factors, however, in that work the reactant's subsystem mode was assumed to be in the low-temperature limit. Although it only makes a minor difference, here, we include as many reactant states as is necessary to converge the rate, and perform the Franck-Condon overlap integrals numerically. Where necessary, we take the continuum states of the Morse oscillator into account, by extending the MLJ formula given in Eq. (15) in the same way as explained for the exact quantum rate [Eq. (A1)] in Appendix A. The reaction rates for this model system computed with various methods as a function of the driving force ε are presented in Fig. 3. The exact (Fermi's golden rule) and MLJ rate calculations are based on the knowledge of the analytic expressions for the energy levels and wavefunctions of the Morse oscillator (see Appendix A). Thus, the only approximation made by MLJ is to treat the bath classically. It can be seen from Fig. 3 that instanton theory, which does not require knowledge of the eigenstates nor even global knowledge of the potential along the subsystem mode, is virtually identical to the exact result when employing a quantum bath with the effective action from Eq. (10). This excellent agreement was expected from the results and analysis seen in previous instanton studies of electron transfer. 70,79,84 When using a classical bath with the action given by Eq. (11), it is slightly less accurate, although then very similar to MLJ theory as they both suffer from the assumption of a classical bath.
The classical golden-rule transition-state theory (TST) rate, outlined in Appendix A, constitutes the classical limit of the quantum rate and would reduce to Marcus theory in the case of a subsystem consisting of displaced harmonic oscillators. The deviation of the TST rate from the exact, MLJ and instanton rates underlines the importance of nuclear quantum effects, which causes the classical rate to differ from the exact results by more than seven orders of magnitude in some cases. The differences are most extreme for the largest driving forces in the inverted regime.
For ε > Λ, the instanton analysis predicts a negative value of τ , which is a clear indication that the inverted regime has been reached, and it requires a subtly different ring-polymer optimization scheme. 84 Despite this, it is noteworthy that the observed turnover in the rate (i.e. the point at which the rates start to decrease with growing driving force) actually occurs at a slightly smaller driving force. This is predicted correctly by all methods tested apart from classical golden-rule TST. In instanton theory, this effect is caused by the prefactor in Eq. (30) as the effective reduced action in the exponential has its minimum at ε = Λ.
The fact that in the inverted regime the MLJ, instanton and exact quantum rates almost coincide, reveals that practically all the quantum effects in this regime originate from the subsystem and not from the bath. Although the quantum subsystem still plays a dominant role in the normal regime, it is clear that there is also a small quantum effect from the bath, which explains the source of the error of MLJ and likewise of rSCI theory when employing a classical bath.
In Sec. V, we will further elaborate on the relationship between rSCI and MLJ theory that is apparent from the results in this section.

C. Comparison with alternative approximations
In this subsection, we will compare the instanton approach with two other approximate methods for including quantum effects into electron-transfer rates, namely the cumulant expansion and the semiclassical Franck-Condon sum. As well as discussing the accuracy of these various methods, we will also focus on the computational effort required for their calculation.
In Fig. 3, we present the rates obtained with the second-order cumulant expansion 118-120 described in Appendix B. To enable a direct comparison with MLJ theory, we employed a classical bath in its calculation, although like with rSCI it would also be possible to use the effective action of a quantum bath. This method is not only commonly used for the study of electron-transfer reactions in anharmonic systems, [121][122][123][124][125][126][127][128] but also for the simulation of optical spectroscopy, 129 vibrational lineshapes 95 and the description of energytransfer processes. 130 In practice often further approximations are invoked to obtain analytical expressions for the rate 131 before the method can be applied to complex problems.
The advantage of the cumulant expansion over an exact (FGR) or MLJ calculation is that it does not require knowledge about the excited state's vibrational eigenstates, but only about its potential-energy surface. However, although the method is exact for displaced harmonic potentials, 95 the results in Fig. 3 clearly demonstrate, that, as opposed to instanton theory, the rates obtained by the second-order cumulant expansion can differ significantly from the exact rates for the Morse oscillator model, with the worst case being at zero driving force in the normal regime. 132 Moreover the rate expression of the cumulant expansion does not satisfy the detailed balance relation for thermal rates, 107,133 in anharmonic subsystems or even in a subsystem of two displaced harmonic oscillators of different frequency.
Here, k 0→1 and k 1→0 are the rate constants of the forward and backward reactions. Note that this relation would normally be written with total partition functions, but here we have already used the fact that in our case Detailed balance is however obeyed by Fermi's golden rule, MLJ theory, all forms of instanton theory and even classical golden-rule TST.
Another method that does not obey detailed balance for anharmonic subsystems is the "semiclassical Franck-Condon sum" (SFC). In fact, the rates computed within this approximation do not even fulfil detailed balance for a subsystem of two displaced harmonic oscillators of the same frequency if the driving force is different from zero. Originally the method was developed to describe spectral line shapes of solids 134,135 and later used to describe electron-transfer in biological systems. 43,136 The derivation of the method for models with a harmonic bath as considered in this paper is outlined in Appendix C. In this case we treat the bath itself within the SFC approximation as this can be done with a closed-form expression. Like the cumulant expansion, it requires knowledge of the vibrational eigenstates of the reactant but not of the product. In accordance with the findings of Siders and Marcus,11,96 it is accurate near the activationless regime, works fairly well in the inverted regime and gives significant errors of more than four orders of magnitude in the normal regime.
Ultimately, the major intrinsic problem of the MLJ method is that it relies on knowledge of the wavefunctions of the reactant and product states and is therefore practically impossible to apply to complex multidimensional problems. The cumulant expansion and SFC methods only go part of the way to improving this situation as they require wavefunctions only for the reactant state. However, numerical integration over the coordinates would still require both potentials to be evaluated over a large grid. Typically therefore at least the reactant potential is approximated by a low-dimensional harmonic oscillator, which introduces an unknown additional error into the predicted rate.
In contrast, instanton optimizations only require information along the tunnelling pathway, which is located close to the hopping point and thus minimizes the computational effort. This advantage of instanton theory over the wavefunction-based methods increases in significance with growing dimensionality of the subsystem. The reason for this is that the instanton pathway always remains one-dimensional, whereas the number of points needed to evaluate the potential-energy surfaces on a grid grows exponentially with subsystem size. It can thus be applied in principle to complex systems without making extra approximations.
It is therefore worth noting that although our instanton approach as well as the SFC method and a number of other theories are labeled "semiclassical", they clearly employ quite different approximations. Not only is semiclassical instanton theory superior in accuracy, it is also applicable to more complex multidimensional anharmonic problems.

V. INSTANTON FORMULATION OF MLJ THEORY
Although MLJ and reduced instanton theory can both be derived from Eq. (5), the resulting methods and rate formulas [Eq. (15) and Eq. (30)] look rather distinct from each other and thus lead to quite different interpretations of the reaction. Marcus-Levich-Jortner theory relies on the wavefunction picture of quantum mechanics and computes the rate as a sum over reactant and product states which will be dominated by one particular reaction channel as shown in Fig. 1(d). Instanton theory, on the other hand, is based on the path-integral formalism of quantum mechanics and is dominated by a path which describes the mechanism during the electron-transfer event.
Another fundamental difference between MLJ theory and the rSCI approach presented in Sec. IV is that, in rSCI, the focus is shifted from the bath modes to the subsystem. The standard MLJ picture as shown in Fig. 1 interprets the reaction in terms of the activation energy in the bath and includes the effect of the subsystem through the shift that they give to the bath potentials. The computation of the reduced instanton approach, however, is carried out directly in the subsystem modes under the influence of the bath. This reflects more appropriately the computational effort put into the calculation of subsystem and bath, as typically the subsystem will be treated in much more detail or on a higher level of theory.
Both interpretations can be useful, but it is not immediately obvious that they can be reconciled, although the common foundation in Eq. (5) suggests that both methods must be related. This idea is reinforced by the fact that the rates obtained for the double Morse oscillator model, shown in Fig. 3, are practically identical when both methods treat the bath classically. In the following we will show that a different derivation of the semiclas-sical instanton approximation leads to an equivalent formulation but which can be used to give the same insights as MLJ theory.

A. Formalism
The objective of this section is to derive an instanton formulation of MLJ theory. The bath is thus assumed to be classical and for simplicity both subsystem and bath are kept one-dimensional here. The formulas do, however, generalize straightforwardly to the multidimensional case.
In order to show the relation with MLJ theory more closely, the convolution formula [Eq. (6)] will again serve as the starting point. The expression for the lineshape function of the bath in Eq. (7b), will be evaluated by a classical phase-space integral, which is one dimensional in both the position and momentum coordinate. After carrying out the integrals in momentum and time, this results in the one-dimensional classical configuration-space integral (40) where the Hamiltonians of the bath [Eq. (3b)] have been replaced by their classical analogues and the independence with respect to τ appears naturally. In addition, we define the potential energy difference in the bath , although we suppress the Qdependence to avoid clutter. For the harmonic bath potential, ∆V b = −2M Ω 2 ζQ = −Λ b Q/ζ. The Q-integral could of course easily be carried out immediately to give the Marcus theory lineshape. However, In order to obtain a picture of the reaction from the point of view of the bath, we leave it for later.
Using this classical result for the bath lineshape function in Eq. (6) and performing the convolution integral leads the approximate rate formula where we define the subsystem lineshape function weighted by the bath thermal distributioñ Note that the effect of the convolution manifests itself in a change of the argument of the subsystem lineshape function I s , which now implicitly depends on the bath coordinate Q via ∆V b . Viewing the expression for the reaction rate with an implicit dependence on the bath coordinates is also the idea that enables the illustration of the MLJ rate by shifted potentials along the bath modes, as shown in Fig. 1. In fact, if Eq. (14) is used for the subsystem lineshape function and the remaining Q-integral over the delta-function in Eq. (40) is taken, the standard MLJ rate formula [Eq. (15)] is recovered.
Here, we seek to treat the subsystem part with semiclassical instanton theory. Note that both the MLJ and instanton version of the subsystem lineshape function emerge from Eq. (8). The difference is induced by the order in which the sums and integrals in Eq. (8) are taken. Whereas in MLJ theory the configuration-space integrals are taken before the sums over states are carried out, in instanton theory these steps are taken in reversed order leading to a path-integral instead of a wavefunction formulation of the reaction rate. Only in the path-integral formulation is it possible to take the steepest-descent integration which leads to semiclassical instanton theory. The instanton subsystem lineshape function is thus given by 137 where again all quantities are evaluated at the stationary point of the exponent S(τ )/ + (∆V b − ε)τ / in the subsystem coordinates q , q and imaginary time τ simultaneously. Using this approximation in Eqs. (42) and (41) defines the instanton formulation of MLJ theory.
Here we show that this approach gives the same result as the reduced instanton theory derived in Sec. IV A. By employing Eq. (42) for the subsystem lineshape function in Eq. (41), the effective action in the exponent becomes Due to the harmonic nature of the bath, the stationary point in the bath coordinates can be solved for analytically. This defines the hopping point at which the electron transfer dominantly takes place. Within the classical limit, it is given by Evaluating Eq. (44) at this point therefore leads S(Q ‡ , τ ) = S r (τ ). So the exponent becomes identical to that of reduced instanton theory (with a classical bath) and hence the value of τ at the stationary point is the same too. The rate expression for this instanton version of MLJ theory is obtained by performing the remaining Qintegral by steepest-descent and using the classical partition function Z b 0 = (β Ω) −1 to give where all quantities are evaluated at the stationary point Q = Q ‡ including the system lineshape function, which implicitly depends on Q through ∆V b . In order to verify the equivalence of this rate expression with Eq. (30), we make use of the rules of consecutive steepest-descent integrations 78,138 Because the spatial subsystem and bath coordinates are independent, the partial derivatives involving Q can be easily evaluated. After rearranging, this results in where ∂ 2 S r ∂τ 2 = ∂ 2 S ∂τ 2 − 2Λ/β . Using these expressions in Eq. (46) shows that this approach is therefore identical to rSCI [Eq. (30)], which is not surprising as all we have done is carry out the same steepest-descent integrations but in a different order.
Following this procedure for the displaced harmonicoscillator models defined in Table I, we obtain the instantons depicted in Fig. 4. Panels (a) and (b) show I(ε − ∆V b ) computed with the semiclassical instanton approximation as a function of Q for the inverted and normal-regime model. As expected, the function is centered around Q ‡ and is well approximated by a Gaussian. From instanton theory, we have therefore obtained a reduced picture of the reaction, but this time the focus is along the solvent coordinate and hence can provide a similar interpretation to that from MLJ theory.

B. Analysis and Mechanistic Insights
In addition to the formal connection between SCI and MLJ discussed in Sec. V A, we will show that, as well as the insight into the tunnelling pathway, it is possible to use instanton theory to extract very similar information about the reaction as is offered by MLJ theory, such as the bath activation energy and the dominant reactant and product vibrational states. We thus suggest that instanton theory may be used instead of MLJ theory for understanding and interpreting electron-transfer reactions in complex anharmonic systems.
In Table II, we present numerical values of the reaction rates for the two models in the normal and inverted regimes models defined in Table I computed with different methods as well as a number of values obtained from the instanton calculation which we will describe later. A comparison of the accuracy of the approaches has already been carried out in Sec. IV B and thus here we simply note a couple of points which are special to this case. The fact that for the inverted-regime model the rSCI rate is even slightly closer to the quantum rate than the MLJ rate can be attributed to a fortuitous error cancellation, as MLJ theory is, in principle, the more accurate method in this case. Because both models consist of displaced symmetric harmonic oscillators, the second-order cumulant expansion is exact in these cases and therefore not shown. In contrast, the rate obtained with the SFC method, while showing decent agreement with the exact rate for the inverted-regime model, exhibits an error of almost one order of magnitude for the normal-regime model. This is in agreement with the findings in Refs. 11 and 96. Just as in the reduced instanton formalism derived in Sec. IV, the instantons computed in the subsystem co-ordinate space, shown in the bottom panels of Fig. 4, consist of paths q n (u) whose energies are not equal but differ by the amount ∆E s . We will show that this energy jump is a good approximation to the difference in energies between the dominant reactant and product vibrational states in the MLJ sum, i.e. E ν 1 − E µ 0 . As explained in Ref. 84, in the normal regime, the trajectories travel in opposite directions away from the hopping point q ‡ , but in the same direction when in the inverted regime. This occurs because τ < 0 in the inverted regime such that the product trajectory travels in negative imaginary time and  Table I. The reduced instanton rates with classical bath and the corresponding values τ at the stationary point of the reduced action S r were obtained from ring-polymer instanton optimizations with 256 beads equally distributed between both electronic states. The contributions to the total effective action from subsystem S n and bath Φ cl are also given. The exact rate is obtained from integration of the flux correlation function of the full system. 78 As described in Appendix C, the SFC approximation is used for both subsystem and bath in order to compute the corresponding rates. For both models, the MLJ rate includes contributions from excited reactant states. All rates, including the Marcus rate for the full system k MT (ε), are given relative to the Marcus rate for the bath of the respective model only k b MT (ε).
which is seen to be equal to ∂ Φ cl ∂τ and just justifies identifying this term as the energy jump in the bath in Sec. IV A. At the stationary point we have ∆E s + ∆E b − ε = 0, which confirms that the total energy is conserved.
As well as predicting the energy jump, we can also predict the reactant and product vibrational states which dominate the MLJ sum. In instanton theory the energies of the two trajectories q n (u) making up the reduced instanton, defined by Eq. (28), indicate the energies with the largest contributions to the thermal rate. In this harmonic system we can relate the energies directly to the vibrational quantum numbers as the energy levels are known. This would of course not be possible in a complex system, although knowledge of the energy in the subsystem before and after the reaction, which provides similar insight, would still be available.
As one can read from the table, for the inverted-regime model, the instanton energies correspond to a transition from the reactant ground vibrational state µ ≈ 0 to the product state ν ≈ 15. The dominant vibrational channel in the normal-regime model is predicted to involve an excited reactant vibrational state µ ≈ 3 and the product state ν ≈ 11. A comparison of these values with Figs. 1(c) and 1(d) reveals that, for both models, the instanton energy picks out the same dominant vibrational channel as MLJ theory. Note that instanton theory does not actually quantize the reactant and product wells as it relies solely on imaginary-time trajectories which exist only in the classically forbidden regions. It does not therefore give integer values for the dominant states. This is however not a serious concern as there is no particular relevance of the individual state with the largest contribution because typically MLJ theory predicts that a cluster of states are involved and thus any prediction within the cluster is practically as good. 139 Furthermore, for a subsystem in conjunction with a classical harmonic bath, the activation energy of the bath modes can be easily recovered using Eq. (45) to give which should be evaluated at the stationary value of τ . The values for the bath activation energy obtained from rSCI theory are also given in Table II and are in good agreement (i.e. with an error less than the thermal energy) with the results obtained from MLJ in theory given in Sec. III B.
In the inverted-regime model, as can be seen in Fig. 4(c), the bath activation energy is thus substantially lower than it would be if there were no subsystem, for which it would correspond to the point where the potentials cross. The presence of the subsystem therefore leads to a significant speed-up of the reaction, which explains why the rate for the full system in Table II is many orders of magnitude larger than the corresponding reaction taking place in the bath only. However, the rate is not only dependent on the bath activation energy but also depends on the action of the subsystem instanton, as can be seen from Eq. (44). As previously discussed, the stationary value of the bath configuration, Q ‡ , is associated with an energy jump ∆V b (Q ‡ )−ε, that must be compensated by ∆E s = ∆V s (q ‡ ) with an equal magnitude but opposite sign in the subsystem in order to satisfy energy conservation. Panels (c) and (e) of Fig. 4 illustrate that minimizing the bath activation energy causes the tunnelling pathway in the subsystem to lengthen which increases the system action. The elongation of the path can be understood from Fig. 4(e) which shows that, in order to reach a point where the potential-energy difference between the subsystem potentials exactly compensates the energy jump in the bath, the system has to travel "uphill" to the left. Hence, in general a compromise has to be made between minimizing the bath activation energy and the subsystem action. Only in one particular case is the magnitude of the energy jump at the reactant mini-mum in the bath (including the driving force) identical to the potential-energy difference at the reactant minimum in the subsystem such that an activationless reaction becomes possible. In this special case, the energy jump in the bath is Λ b − ε and in the system is ∆V s (q (0) min ), which for our example where V s 0 (q (0) min ) = V s 1 (q (1) min ), leads to the requirement ε = Λ. The rates in the inverted regime are much faster than in the fully classical treatment, because the reaction within the subsystem can proceed via quantum tunnelling, as depicted in Figs. 4(e) and (f), instead of relying on thermal activation. Altogether, this implies that, although the turnover curve in the inverted regime is not as steep as it would be according to the classical theory, it does not become independent of ε.
On the other hand, in the normal regime, the bath activation energy in Fig. 4(d) is seen to be higher than it would be without the presence of the subsystem. This causes the rate for the full system to decrease relative to the electron-transfer reaction in the bath only. The tunnelling effect increases the rate relative to a classical calculation, although typically not as dramatically as in the inverted regime. This can also be understood from an analysis of the instanton tunnelling trajectories as was explained in Ref. 84.
Thus, we were able to show how practically all insights from MLJ theory including, first and foremost, the dominantly contributing reactant and product energies can equally be obtained from reduced instanton theory. Semiclassical instanton theory further allows one to attain this understanding of the reaction even in complex, anharmonic systems. This information is complemented by the localization of the optimal tunnelling pathway in the subsystem, which can be interpreted as the reaction mechanism in configuration space.

VI. CONCLUSIONS
We have developed an instanton formulation of MLJ theory, which focuses on the subsystem while including a classical or quantum harmonic bath implicitly. This provides a practical method to complement the simulation of electron-transfer reactions of multidimensional anharmonic subsystems by the effect of a solvent environment. Thus, the method is ideally suited to study problems that necessitate multiscale modelling.
Electron-transfer rates have been calculated and compared to results from several other methods for an asymmetric anharmonic model and the results demonstrate that reduced semiclassical instanton theory is in excellent agreement with either the exact rate or the MLJ rate depending on whether the bath is assumed to be classical or not. Thus we argue that semiclassical instanton theory can be reliably employed in situations which have previously been simulated by MLJ theory.
In addition to MLJ theory, we have also compared our approach to the second-order cumulant expansion, a popular method commonly used to describe electron-transfer and optical transition rates, and to the semiclassical Franck-Condon sum. The results obtained with both these approximations exhibit severe errors, especially in the normal regime and in fact, unlike instanton theory, neither the cumulant expansion nor the SFC approximation satisfy the detailed balance relation [Eq. (39)]. This underlines the fact that, although both the SCI and SFC methods have been termed "semiclassical", the approximations are quite unrelated.
We also compared and contrasted the insight that MLJ and instanton theories can offer into the mechanism of electron-transfer reactions. The traditional MLJ picture is shown along the bath coordinates, in which the subsystem has an effect by shifting the reactant and product potential by their respective internal energy levels. Although undoubtedly simple and intuitive in one dimension, this picture quickly becomes convoluted when a multidimensional anharmonic subsystem has to be considered. There is also little insight given into the tunnelling dynamics of the subsystem itself.
Instanton theory, on the other hand, automatically locates a unique reaction coordinate which describes the optimal tunnelling pathway of the subsystem modes. In analogy to the dynamics of open quantum systems, the addition of a bath changes the instanton pathway in the subsystem such that, due to energy exchange between subsystem and bath, the reduced instanton exhibits an additional jump in energy at the hopping point. Although the energy in the subsystem is therefore not conserved by the electron-transfer reaction, the excess energy is absorbed by the bath such that the total energy is conserved as it of course should be. This picture of tunnelling under the barrier along a reaction coordinate reflects the typical situation of practical simulations, where the focus is on the subsystem under the influence of a surrounding solvent bath.
Nonetheless, we have also discussed how instanton theory is connected to MLJ theory by deriving them both from a common expression. This shows that in principle similar insights can be extracted from either method. In particular, we show that instanton theory can successfully predict the same dominant initial and final vibrational state of the system before and after the electrontransfer event as MLJ theory.
Instanton theory overcomes the main disadvantage of MLJ theory, which is that it requires knowledge of the energy levels and wavefunctions of the subsystem. Because of this, applications of MLJ theory are often limited to a harmonic-oscillator approximation, which introduces an uncontrolled error when simulating an anharmonic system. Hence, although rSCI (with a classical bath) is technically an approximation to MLJ theory, in many anharmonic cases it will lead to more accurate results due to its ability to account for anharmonicity along the tunnelling pathway. In conjunction with a ring-polymer discretization, instanton theory can be applied directly to multidimensional anharmonic problems. The application of this theory to electron-transfer reactions, spin transi-tions and energy-transfer processes of molecular systems in combination with high-level ab-initio electronic structure methods, as has been used in previous instanton studies, [92][93][94] will be integral part of future work.

ACKNOWLEDGEMENTS
This work was financially supported by the Swiss National Science Foundation through SNSF Project 175696.
Appendix A: Quantum and classical rate formulas for an anharmonic subsystem mode in conjunction with a harmonic bath In order to put the instanton and MLJ results shown in Fig. 3 into context, we also present the exact quantum rates for this system and their classical limits. This enables us not only to directly check the quality of the results obtained with the approximate methods, but by comparison with the classical rates also allows an estimation of the relevance of nuclear quantum effects. In this section, we harness the formal framework laid out in Sec. II to derive the required rate formulas making use of the fact that we can analytically integrate out the coordinate-dependence of the harmonic bath.
If the wavefunctions of the subsystem are known, as is the case for the two crossing Morse oscillators used in Sec. IV B, the trace over the subsystem degrees of freedom in Eq. (7a) can be evaluated exactly in the wavefunction representation. Thus, the exact quantummechanical rate can be computed by the formula where sums are taken over the bound states of the Morse potential and integrals are carried out over the energies of the energy-normalized continuum states. Expressions for the wavefunctions can be found in Refs. 140 and 141. Numerical integration was used to obtain the Franck-Condon overlaps and to perform the integral over time. In order to make the latter converge easily, the imaginary-time variable τ was chosen appropriately (i.e. using the value obtained from the instanton optimization).
The classical limit of this rate can be obtained in a similar way except that the trace in Eq. (7a) is evaluated by a classical phase-space integral 106 and the classical limit of the effective bath action is used [Eq. (11)]. For a one-dimensional subsystem, this gives the classical golden-rule transition-state theory rate where the reactant partition function is computed by a classical phase-space integral. The remaining integral in the subsystem mode can either be taken numerically, as was done to generate the results in Fig. 3, or by steepest descent, which would be an excellent approximation in this case. The quantum-mechanical rate [Eq. (A1)], as well as the MLJ rate [Eq. (15)] correctly reduce to the TST expression in Eq. (A2) in the high-temperature or low-frequency limit, while instanton theory [Eq. (30)] reduces to the steepest-descent version of it.
In the special case that the subsystem consists of displaced harmonic oscillators, Eq. (A2) reduces to the Marcus theory expression [Eq. (13)] except that in this case the reorganization energy should be the sum of the subsystem and bath reorganization energies.

Appendix B: Cumulant expansion
Another approximate way of computing correlation functions and therefore also to calculate electrontransfer reaction rates is the so called "cumulant expansion", 118-120 which in our formulation will be applied to the lineshape function of the subsystem Eq. (7a).
In its conventional formulation τ is set to zero and we rewrite Eq. (7a) as where the correlation function is R(t) = (Z s 0 ) −1 Tr s e −(β −it)Ĥ s 0 / e −itĤ s 1 / .
The time-dependent terms inside the trace can equally be rewritten as a time-ordered exponential according to e +itĤ s 0 / e −itĤ s 1 / =T e −i t 0 ∆V s I (t )dt / whereT is the time-ordering operator and we make use of the interaction picture to give ∆V s I (t) = e +iĤ s 0 t/ ∆V s e −iĤ s 0 t/ , where ∆V s =Ĥ s 1 −Ĥ s 0 = V s 1 (q) − V s 0 (q). 95 This exact expression can then be expanded in a time-ordered power series with respect to ∆V s I . Motivated by the analytic solution for the correlation function of a system of displaced harmonic oscillators [Eq. (9)], one makes the ansatz R(t) = exp [−Γ(t)] where the exponent is defined as a sum of cumulants where Γ j (t) is of jth order in ∆V s I . Comparing the two expansions, the first two terms in Eq. (B3) are given by 95 where ∆V s µµ = ∞ −∞ ψ µ 0 (q) * ∆V s (q)ψ µ 0 (q) dq and these integrals over the one-dimensional subsystem coordinate are evaluated numerically. The terms in the sum of Eq. (B5b) with µ = µ can be evaluated by L'Hpital's rule. For the case of displaced harmonic oscillators, this expansion of the correlation function up to second order gives the exact result, as all higher order terms vanish. 95 In the general, anharmonic case, however, the quality of the approximation is unclear. One could of course extend the method to higher orders, but the series is unlikely to converge quickly to correct result.
The computational advantage of the cumulant expansion over the golden-rule formula is that only the eigenstates of the reactant electronic state need be known. It is thus perhaps most useful when computing absorption spectra from a ground electronic state to an excited state, for which the ground state is well approximated by a harmonic oscillator, but not the excited state. However, a significant knowledge of the product potential-energy surface is still required in the region where the wavefunction overlaps in Eqs. (B5) are sizeable, which can be expensive to compute.
This result for the subsystem's lineshape function can be easily combined with the lineshape function of the harmonic bath from Eq. (9) by performing the convolution integral in v and integrating over the resulting delta function. The rate expression based on the second-order cumulant expansion for the subsystem part is therefore given by where either a quantum or classical bath can be employed by using the respective expressions for the actions in Eqs. (10) and (11) and the final time-integral is carried out numerically. Note that this cumulant expansion leads to a completely different approximation from that of Wolynes theory 59 even though the latter can also be thought of as a type of cumulant expansion. In contrast to the approach described here, Wolynes theory carries out the time integral by the method of steepest-descent and computes the short-time limit of the correlation function by path-integral sampling. For the systems studied in this work, Wolynes theory would give similar results to those of instanton theory (identical in the case of a harmonic system), although for certain more complex systems it has been shown to break down. 70,85 Unlike the cumulant expansion and instanton theory, 84 it is also not directly applicable to the inverted regime, although an extrapolation method which extends it in this way has been suggested. 67 D j=1 Λ b j γ j coth γ j /Λ b , which is defined in terms of the reorganization energy associated with a single bath mode Λ b j = 2M Ω 2 j ζ 2 j and γ j = β Ω j /2.
Following the formalism laid out in Sec. II and performing the convolution integral in v first, this leads the rate equation where the integrals over the anharmonic subsystem mode have to be carried out numerically.
In the special case in which all modes are displaced harmonic oscillators, all degrees of freedom can be assigned to the bath. Then, the rate formula is directly given by k SFC (ε) = ∆ 2 2 I b SFC (ε). It is easy to see that this is in error because it predicts results symmetric around ε = Λ, whereas the true result is known to be significantly skewed unless in the classical limit. 11,20,51 One way to understand the causes of this error has been explained in terms of WKB theory in Ref. 142.