Revisting Nuclear Tunnelling in the Aqueous ferrous-ferric Electron Transfer

The aqueous ferrous-ferric system provides a classic example of an electron-transfer processes in solution. There has been a long standing argument spanning more than three decades around the importance of nuclear tunnelling in this system, with estimates based on Wolynes theory suggesting a tunnelling factor of 65, while estimates based on a related spin-boson model suggest a smaller factor of 10-36. Recently, we have shown that Wolynes theory can break down for systems with multiple transition states leading to an overestimation of the rate, and we suggest that a liquid system such as the one investigated here may be particularly prone to this. We re-investigate this old yet interesting system with the first application of the recently developed golden-rule quantum transition state theory (GR-QTST). We find that GR-QTST can be applied to this complex system without apparent difficulties and that it gives a prediction for the tunnelling enhancement factor 6 times smaller than that from Wolynes theory. The fact that these theories give different results suggests that although it is well known that the system can be treated using linear response and therefore resembles a spin-boson model in the classical limit, this approximation is questionable in the quantum case. It also intriguingly suggests the hypothesis that the previous predictions were overestimating the rate due to a break down of Wolynes theory.


Introduction
The realms of chemistry and biology serve us with a colourful variety of reactions affected by nuclear tunnelling 1 . In chemistry, tunnelling is predicted to be important in reactions under a wide range of conditions from astrochemical reactions occurring on cosmic dust 2 and nuclear fusion in stars 3 to organic chemistry, where even heavy-atom tunnelling has been identified 4,5 . Biological systems have also been suspected of employing nuclear tunnelling, for instance in photosynthesis taking place in bacteria 6 or during enzyme catalysis 7,8 . To resolve such controversial hypotheses, a reliable method to calculate effects of nuclear tunnelling is clearly desirable. Such a theory will be useful to quantify the relevance of tunnelling in a given reaction.
In this work we focus specifically on the example of electrontransfer reactions [9][10][11] . These reactions are nonadiabatic and governed by a change of electronic state and therefore cannot employ the Born-Oppenheimer approximation. 12 The rate is however well described by Fermi's golden rule 13,14 , although in practice this cannot be evaluated for complex molecular systems as it requires complete knowledge of the internal eigenstates of the a Laboratory of Physical Chemistry, ETH Zürich, 8093 Zürich, Switzerland † These authors contributed equally * E-mail: jeremy.richardson@phys.chem.ethz.ch system.
Modern quantum rate theories 15 are typically based on the path-integral approach to quantum mechanics 16 , which allows tunnelling and other nuclear quantum effects to be included efficiently into molecular simulations using a quantum-classical correspondence 17 . However, because the rate is not defined as a simple expectation value of the density matrix, but rather in terms of a time correlation function, 18 it is by no means trivial to calculate rates in this way and further approximations are required. In this paper, we will concentrate on quantum transition-state theories and not consider dynamical methods 19,20 .
Semiclassical instanton rate theory [21][22][23] predicts the tunnelling rate and mechanism via locating the optimal tunnelling pathway (the instanton) based on a stationary-action principle. Based on a similar first-principles derivation as in the normal regime 24,25 , instanton theory has been extended to treat electron-transfer reactions [26][27][28] in both the normal and inverted regimes. 29 It has the most rigorous derivation of the methods discussed in this paper, shows excellent agreement with exact methods on model systems and is well suited for gas-phase electron-transfer reactions. However, for liquid systems, it is formally not valid to apply instanton theory 25 , although in some cases approximate application is possible by using an implicit solvent model or by freezing all atoms not expected to be involved in tunnelling at the transition state (TS) geometry 30 . For the general case, an extension of instanton theory that allows for sampling is desired.
Wolynes theory 31 is an approximate quantum rate theory which describes electron transfer in Fermi's golden-rule limit. It is defined in terms of path integrals which can be evaluated using an N-bead discretization with each bead assigned to either the reactant or product electronic state. The method of path-integral molecular dynamics (PIMD) 32 opens Wolynes theory up to the sampling tools of molecular dynamics (MD) calculations and accordingly makes it a computationally feasible theory for atomistic systems. [33][34][35] . Lawrence and Manolopoulos have recently shown that Wolynes theory can also be successfully extrapolated to the Marcus inverted regime 36 . Wolynes theory has been thoroughly investigated not only for atomistic but also for model systems such as the spin-boson model, where it compares very well to exact results 33,37 , because it recovers the stationary phase-approximation 38 . A limitation to Wolynes theory however is, that it does not tend to the classical limit for anharmonic systems 27,39 . Recently, we have also pointed out another crucial limitation of Wolynes theory, which is that its approximations break down when a system consists of two or more transition states. 40 This break-down can manifest itself as an overestimation of the reaction rate by more than an order of magnitude in either the classical or the quantum limit. This may lead to the prediction of an artificial tunnelling factor. The break-down of Wolynes theory can be related to its lack of connection to instanton theory, as it is observed that paths far from the instanton are sampled. In systems where this break-down is observed, Wolynes theory does not necessarily sample paths close to the diabatic crossing seam but can rather include unphysical configurations far from the seam. This makes any mechanistic insight or a correct rate prediction impossible.
The quantum-instanton method 41 suffers in a similar way when applied to strongly asymmetric barriers, which can be explained from an analysis in terms of semiclassical pathways and corrected by introducing a projection to connect it to the instanton. 42 A further example to back this line of argumentation is the success of ring-polymer molecular dynamics (RPMD) [43][44][45] , which was shown to be closely connected to the semiclassical instanton rate theory in the deep-tunnelling regime 46 . Building on RPMD rate theory in the adiabatic limit [47][48][49] , attempts were made to extend it to treat the nonadiabatic limit. 15 Two such attempts are the kinetically-constrained RPMD 50 and the isomorphic RPMD method, which do not always give reliable tunnelling factors. 51 One can in turn relate this behaviour to their lack of connection to instanton theory.
We therefore proposed golden-rule quantum transition-state theory (GR-QTST) 39 in order to overcome issues of possible break-down behaviour by keeping a relation to instanton theory, but also retain the advantageous feature of Wolynes theory which includes not only the instanton but also paths in its vicinity. This method is computed in a similar way to Wolynes theory, except that a constraint is imposed on the sampled paths such that the energy on the reactant and product states must match. Adding such a constraint has been proposed as a general approach for defining quantum transition-state theories. 52 This constraint is automatically obeyed by all instantons, which ensures a strong connection to instanton theory, and it also retains the correct classical limit. GR-QTST has been shown to perform very well for model systems in both the normal and inverted regimes, 39 including the multidimensional spin-boson model. GR-QTST was also investigated for systems with multiple transition states, where Wolynes theory breaks, and provides accurate rate predictions 40 . For the systems tested so far, we claimed it was the most accurate imaginary-time path-integral method currently available. However, GR-QTST has not previously been applied to atomistic simulations. This work therefore aims to investigate the applicability of GR-QTST as well as to see what physical insights it can offer by revisiting the early papers on the aqueous ferric-ferrous electron transfer 33,53,54 . This is a prototypical atomistic system for which a computationally inexpensive force field is readily available 33,53 . The seemingly simple interactions in this system forge a rough, high-dimensional, anharmonic potential energy surface (PES), and display high levels of complexity due to it being atomistic, which is more realistic and complex far beyond any of the models previously studied by GR-QTST.
The system is depicted in Fig. 1, despite its seemingly innocent appearance, there has been a long standing argument over the importance of nuclear tunnelling in this system at room temperature. The nuclear tunnelling enhancement factor reported by Bader et al. 33 and calculated based on similar ideas to Wolynes theory is approximately 65, suggesting a significant contribution from nuclear tunnelling. This estimate is significantly larger (∼ 6 times) than other predictions made at the time 55 . The controversy over the correct rate prediction for this system is interesting in itself, but due to the range of suggested tunnelling enhancement factors the aqueous ferrous-ferric system is a good atomistic test case especially to investigate in view of possible overestimation of nuclear tunnelling enhancement and the applicability of Wolynes theory. Therefore, naturally the aqueous ferric-ferrous system is worth revisiting with newly developed rate theories.
In order to reexamine the earlier findings of Chandler and coworkers 33,53 and add the investigation of GR-QTST for this system, we recapitulate the various rate theories under study in Section 2. The computational details of the implementation of each rate theory are given in Section 3 and the results are given and discussed in Section 4. We conclude in Section 5 on the quality and appropriateness of the various quantum rate theories and Fig. 2 Illustration of a ring polymer on two diabatic PESs V 0 and V 1 in a two-dimensional nuclear configurational space. Only contours for the lowest PES are shown at any configuration. The blue (red) beads of the ring polymer represent the imaginary-time path on the reactant (product) electronic state. In this example, N 0 = 6 and N 1 = 4 giving a total of N = 10 beads. discuss our work in the context of that of others.

Theory
The quantum Hamiltonian describing an electron-transfer problem is 12Ĥ in which ∆ is the electronic coupling, andĤ n = ∑ D j=1p 2 j /2m j + V n (x) is the nuclear Hamiltonian for the electronic state |n with the PES V n (x), where x = (x 1 , ..., x D ) is the nuclear geometry and the index j runs over each of the D nuclear degrees of freedom of the system with momentum p j and associated mass m j . Following the work of Bader et al. 33 , ∆ is assumed to be a constant, which is known as the Condon approximation. The rate, in the limit of small ∆, is in principle given by Fermi's golden rule 13 , which is commonly approximated using Marcus theory 10 . For a symmetric system, the rate is given by a simple equation: where Λ is the reorganisation energy, defined by the average energy gap between the two PESs for a classical ensemble in the reactant state. Despite its successful use in a wide range of applications, it ignores nuclear quantum effects. Various quantum methods for evaluating electron-transfer rates in the golden-rule limit have been derived, 14,56,57 but in this paper, we shall only consider methods based on the imaginary-time path-integral formulation, which are applicable for complex molecular systems described by atomistic Hamiltonians. In electron-transfer theory, cyclic paths are formed from joining together an open-ended path on the reactant state and an open-ended path on the product state. For a ring-polymer representation of these paths, we introduce λ as a dimensionless order parameter, which determines how the ring-polymer beads are distributed between the two diabatic states |0 and |1 . It is defined by where N 0 , N 1 and N are integers according to the discrete distribution of beads, i.e. N 0 on the diabatic state |0 and N 1 on the diabatic state |1 . The total number of ring-polymer beads is N.
An illustration of such a ring polymer is given in Fig. 2. The two extreme distributions assign all beads to just one diabatic state can be described with the order parameter λ as λ = 0 for the case that all beads are on the reactant PES V 0 and λ = 1 for the case that all beads are on the product PES V 1 . The unconstrained ensemble of ring polymers can be sampled using thermostatted PIMD based on the following Hamiltonian: where ω N = 1/β Nh with β N = β /N and β = 1/k B T , k B is the Boltzmann constant, T is the temperature and x = {x (1) , . . . , x (N) } are the positions of the beads with conjugate momenta p. The cyclic index i runs over each bead such that x (0) ≡ x (N) . Unconstrained free energies are defined in terms of the ensemble of ring polymers by RP dx dp (5) and the reactant free energy by F 0 = F u (0).
First we define Wolynes rate theory 58 , which was derived via a second-order cumulant approximation of the time-correlation function 59 and is given by where F u (λ * ) is the maximum unconstrained free energy with respect to the order parameter λ . For symmetric systems, the maximum occurs at λ * = 0.5.
The GR-QTST method approximates the golden-rule rate using the ansatz 39 : in which λ * is the same as that used in Wolynes theory, i.e. the maximum of the unconstrained free energy according to the approach introduced by us in Ref. 40. F c (λ ) is the free energy under the constraint σ λ (x) = 0 given by The virial energy estimators of the product and reactant paths, E v 0 and E v 1 , are defined as in Ref. 40 and are functions of the potentials and gradients of the beads corresponding to one particular state. The constraint is designed to enforce energy conservation for the ring polymers sampled in the simulation, which is known to give a strong connection to quantum transition-state theories. 52 The classical rate in the golden-rule limit is defined by 12,52 where F cl 0 is the classical free energy of the reactant, and F cl c is the classical free energy of the system constrained at the crossing seam, defined as where There are thus conceptual differences between Wolynes theory and GR-QTST. Wolynes theory relies on a steepest-decent approximation to the time integral of the flux-flux correlation function in the golden-rule limit 31 , whereas GR-QTST incorporates the physical requirement of energy conservation enforced by the virial energy estimator 39 . Both methods are approximations to the true quantum rate, but can be shown to be very accurate for simple systems such as the spin-boson model. 39 It can also be shown that GR-QTST reduces to the classical rate, Eqn. (9), in the high-temperature limit of any system when the ring polymers collapse. 39 However, the same is not necessarily true of Wolynes theory. 52 In particular, we have shown that Wolynes theory can break down for systems with two or more different transition states, due to the fact that only one λ * value is used which cannot simultaneously be optimal for all transition states. In these cases at least, GR-QTST is expected to be more accurate as its rate is approximately independent of the choice of λ * and employs the energy constraint to ensure the correct sampling of each transition state. 40

Computational Methods
We computed rates from Wolynes theory and GR-QTST in the classical limit according to Eqns. (6) and (7). The free-energy term in the Wolynes rate was calculated using thermodynamic integration along the order parameter λ , where the ensemble average of the free-energy derivative, can be obtained from sampling an unconstrained ring-polymer ensemble with the Hamiltonian H RP and histogramming the probabilities of sampling a specific value of the function σ λ * (x), which are defined by The constrained free energy can then be expressed in terms of the sampling probability as This procedure is computationally feasible if the unconstrained simulation samples enough configurations which obey the constraint (σ λ * (x) = 0). If this condition is not fulfilled, the δ -TI method as described in Ref. 40 could be applied in combination to calculate F c (λ * ). However, this was not necessary for the system studied in this work as can be seen from the histograms in both the classical and the quantum limit shown in Fig. 3, which are peaked around σ = 0. Due to this connection in methodology, one can use unconstrained PIMD simulations as a first step towards either Wolynes theory or GR-QTST rate calculations.
Both the Wolynes and GR-QTST rates can also be evaluated in the classical limit. The rate expressions k cl Wolynes and k cl GR-QTST are very similar to Eqn. (6) and (7), with the only difference being that the employed free energies are replaced by their classical counterparts F cl u , F cl c and F cl 0 . In the case of the classical Wolynes rate, the corresponding free-energy difference, F cl u − F cl 0 , can be computed analogously to the quantum case (Eqn. 11) with the ring polymer collapsed onto a single classical particle defined by the Hamiltonian mass scale μ order parameter λ apply constraint Fig. 4 Thermodynamic cycle containing all the free energy integration schemes of Wolynes theory and GR-QTST in the quantum and classical limit. Constrained free energies are obtained at a value of λ = λ * .
Alternatively, F cl u − F cl 0 can be calculated by scaling up the mass of all the atoms (by multiplying m j by µ and taking the limit µ → ∞), which in effect collapses the ring polymer, making it behave classically. We found that both methods give results within each other's error bars for this system and in order to obtain a freeenergy change we employed the mass scaling method.
Note that k cl Wolynes is not the same as the classical rate expression k cl in Eqn. 9, because Wolynes theory does not necessarily tend to the correct classical limit for general systems. However, when dF cl u dλ is a linear function of λ (i.e. for spin-boson models), one can show that k cl Wolynes is the same as the rate of Marcus theory (Eqn. 2), by plugging in dF cl u dλ = Λ(1 − 2λ ) into Eqn. 6. In contrast, as we have shown in Ref. 39, GR-QTST reduces to the correct classical expression, Eqn. (9), in the high-temperature limit where the ring polymer collapses. This gives us two other possible methods for calculating the classical rate, either by sampling H (λ * ) cl to obtain F cl c using the same procedure given by Eqn. 14, or alternatively by scaling up the mass of all the atoms in a GR-QTST simulation.
In order to understand and to visualise the relation between all the different free energy terms in the rate theories introduced above, we constructed a thermodynamic cycle as shown in Fig. 4.
The free-energy calculations necessary to compute Wolynes theory form the top and bottom horizontal legs of the left side of the thermodynamic cycle and relate the reactant (λ = 0) to the stationary point (λ = λ * = 0.5) in both the quantum and the classical limit. In order to compute GR-QTST, one additionally needs the free-energy calculations corresponding to the horizontal legs on the right side of the thermodynamic cycle. They can be viewed as an extension to the Wolynes free-energy calculations.
Each vertical leg in the thermodynamic cycle represents a thermodynamic integration from the quantum nuclei to the classical limit using the mass scaling factor µ as the order parameter (mass-TI) [60][61][62][63] . This means that each leg of the thermodynamic cycle on the left can be calculated from independent simulations. Therefore, we use the left thermodynamic cycle to validate the accuracy and reliability of the free-energies obtained from classical and quantum simulations, giving us further confidence in the Wolynes rates we computed. The free-energy differences, ∆F 0 and ∆F u were calculated by performing two sets PIMD simula-tions with H (λ ) RP at λ = 0 and λ = λ * . The thermodynamic integrand was obtained for 10 different mass-scale factors from µ = 1 up to µ = 100. The contributions from larger µ values are also accounted for via a coordinate transform in the thermodynamic integration 61,62 . We do not compute ∆F c as we found it to be numerically unstable to calculate due to the fact that the virial kinetic energy estimator is not valid for constrained PIMD simulations. This free-energy change can however be inferred by completing the cycle.
As well as the absolute reaction rates, another quantity of great interest is the tunnelling factor Γ, defined as the ratio of the quantum and the classical rate. It should be pointed out however, despite what the name suggests, Γ includes all nuclear quantum effects and not only tunnelling effects. We define Γ for each quantum rate theory, for example where the exponent is defined by as can be seen from Fig. 4. Defining the tunnelling factor in this way gives the most fair comparison between methods as it allows for error cancellation in the case that Wolynes theory breaks down and overestimates the rate in both the classical and quantum limit and may cancel. A relative tunnelling factor as defined by Γ therefore describes the tunnelling enhancement factor described by a given theory and avoids inconsistencies by cross-comparison of different theories. The prefactor ratio A Wolynes /A cl Wolynes is generally not exactly equal to 1, although the equivalent term for GR-QTST is identically 1 and does not therefore appear. Accordingly, the tunnelling factor for GR-QTST can be defined in the same way as Γ GR-QTST = e −β ∆∆F c with the exponent hence defined as The interactions in the aqueous ferric-ferrous system are defined by the interatomic forces and pseudopotentials described in Ref. 53 with the exception of the water model. In contrast to the formerly used rigid single point charge (SPC) water model 33,53,64 , we apply the flexible q-TIP4P/F water model 65 , which was specifically developed to suit PIMD simulations 65 . In particular, it can correctly capture the delicate balance between the competing quantum effects in water, compared to rigid or harmonic water models 65,66 . The reorganisation energy calculated with this water model is 108.5 ± 0.9 mHartree (68 kcal mol −1 , 2.95 eV). This is fairly close to the result from the previous work 33, 53 using a rigid water model, which found a reorganisation energy of 128 mHartree (80 kcal mol −1 , 3.5 eV). The water model we applied also gives a reorganisation energy slightly closer to the experimental estimate of 2.1 eV 55,56 . We also tested the sensitivity with respect to the oxygen-ion interatomic potential, and found no qualitative impact on our findings even when we reduce that by a factor of two.
The q-TIP4P/F and SPC water models include polarisation ef- fects in a mean-field way 67,68 and hence belong to the class of non-polarisable water models, which are computationally affordable and allow for extensive simulations. The application of an explicitly polarisable water model is crucial to describe effects in surface chemistry and clusters 69,70 . The explicit treatment of polarisation is expected to lower the estimate of the reorganisation energy also in the aqueous ferrous-ferric system [70][71][72] and it is therefore not without controversy to employ a non-polarisable water model. Ultimately, the choice of water model will of course affect the quantitative results, but will not hamper our ability to compare the different quantum rate theories. Nevertheless, it should be pointed out that several suggestions on polarisable water models and improved treatment of the solvent models were made in the literature 55,[67][68][69][70]73 .
The ion-ion distance was treated using a fixed-atom implementation at an interionic distance of r = 5.5 Å, which was determined to be the most probable interatomic distance for electrontransfer reactions 33,53 .
In order to obtain the rates for Wolynes theory and GR-QTST, we performed a set of PIMD simulations with N = 24 ring polymer beads, where each set contains 10 starting configurations. The starting configurations were picked randomly from a long MD simulation of the system with 265 water molecules in a cubic box of box length 20 Å (to give a water density of 10 3 kg m −3 ) using periodic boundary conditions. The temperature was set to 300 K and kept constant using the Andersen thermostat. Each set of ring polymers was then sampled under these conditions for 44,000 steps at 13 values of the order parameter λ (see Eqn. (3)) for each integral N 0 value in the range N 0 ∈ [0, 12] to perform a thermodynamic integration along the order parameter λ . The only additional information required to obtain a GR-QTST rate from such a PIMD simulation is a histogram of the the energy constraint σ λ * (x). This effectively requires no extra cost as the energy constraint calculation σ λ * only includes potentials and forces on the PES and only has to be calculated for the stationary point λ * . In terms of statistical error, it suffices to calculate the energy constraint function for 24,000 steps and 4,000 steps are used for equilibration. The mass-TI calculations were accordingly also sampled for 24,000 steps and 4,000 equilibration steps. The timestep was chosen to be 0.5 fs and the equilibration time was set to 2 ps (4000 steps) in both cases. As we have chosen this setup in close analogy to the setup of Ref. 53, the same considerations in terms of finite size effects and potential cutoffs apply. All of the classical MD forces are calculated using LAMMPS 74 .

Results and Discussion
Our aim is to quantify the quantum effects present in the aqueous ferrous-ferric system and thereby to address the controversy of the magnitude of the effect of quantum tunnelling on the reaction rate. In this section we present results from both Wolynes theory and our newly developed GR-QTST 39,40 and discuss the predictions for rate constants and isotope effects from these two different approaches. We investigate possible pitfalls of each theory and discuss how they affect the rate of this system. We then revisit the earlier studies 33  nelling enhancement factors presented for the aqueous ferrousferric system in the literature 33 . Underlying all the rate calculations are free-energy differences which were defined in the thermodynamic cycle introduced in Section 3. In Fig. 5 the results of these free-energy differences obtained from our simulations are given. Before comparing the results obtained from the two quantum rate theories, we first check for consistency of the calculations of the various free-energy paths presented in Fig. 5. Note that the Wolynes rate is defined in terms of the quantities belonging to the thermodynamic sub-cycle on the left, whereas GR-QTST depends also on those on the right. The free-energy change ∆∆F u as defined in Eqn. (17) can be obtained either as a difference of a thermodynamic integration in the quantum and the classical limit, which amounts to (F cl u − F cl 0 ) − (F u − F 0 ) = 3.9 ± 0.2 mHartree, or the difference between the mass-TI calculations at the reactant and stationary-point ensembles ∆F u − ∆F 0 = 3.6 ± 0.3 mHartree. The free-energy differences ∆F 0 and ∆F u obtained from the different mass-TI calculations are significantly larger in magnitude, because they include the change from classical to quantum nuclei and therefore include the zero-point energy of the system. The consistency (within error bars) of the free-energy difference ∆∆F u calculated via the two different routes also means that the Wolynes sub-cycle is closed, which confirms that our simulations are converged. In order to avoid the numerical errors inherent to a subtraction of large numbers we use the leg of the thermodynamic cycle employing thermodynamic integrations in the classical and quantum limit for further calculations.
Rates are then calculated according to Eqns. (6), (7) and (9) from the changes in free energy and are listed in Table 1 for both the classical limit (µ → ∞) and the quantum limit (µ = 1). It is interesting to note that the quantum rate predictions of GR-QTST and Wolynes theory do not agree. Both methods have been tested Table 1 Calculated rates in atomic units in the classical and quantum limit using different rate theories. The classical Wolynes rate calculated with collapsed ring-polymer is almost identical (7.2 ± 1.0 × 10 −11 ). Wolynes 7.6 ± 1.0 × 10 −11 5.3 ± 0.7 × 10 −9 2.7 ± 0.3 × 10 −9 k GR-QTST 6.3 ± 0.9 × 10 −11 7.6 ± 0.9 × 10 −10 4.5 ± 0.4 × 10 −10 on the spin-boson model and give excellent and practically identical predictions of the quantum rate. 39 However, due to the conceptual difference of the two theories, in more complex systems one cannot generally expect Wolynes theory and GR-QTST to predict similar rates. This therefore implies that the aqueous ferricferrous electron-transfer reaction is fundamentally more complex than the spin-boson model. A second obvious conclusion is that since the two theories do not agree, at least one rate prediction must be inaccurate. In the following we analyse the two methods to discuss these points.

Making or breaking of Wolynes theory
As presented in Ref. 40, a break-down of Wolynes theory occurs when the system under investigation exhibits multiple distinct transition states which can lead to an overprediction of the rate by orders of magnitude. There are a number of criteria that serve as indicators to identify whether Wolynes theory is applicable, although the absence of these features does not exclude the possibility of at least a minor break down of Wolynes theory. The first criterion is whether Wolynes theory tends to the correct classical limit as the masses are scaled up. Second, one should investigate the distribution of σ λ * (x) which is sampled. A multi-modal distribution would suggest that Wolynes theory is not applicable. Likewise, the distribution of instantons can be used to obtain a qualitative picture on the nature of the transition states in the system and the appropriateness of Wolynes theory.
In the classical (high-temperature or heavy mass) limit, in contrast to Wolynes theory, GR-QTST is known to tend to the correct classical rate 39 , and our simulations are in agreement with this (see Table 1). As suggested above, there is no rigorous argument that requires Wolynes theory to correctly predict the true classical rate and this thus provides a good check that the method gives physically sensible results. In this case, it gives the correct result within the error bars. The success of Wolynes theory in the classical limit can be related to the approximately linear behaviour of the free-energy derivative with respect to the order parameter, λ , as shown in Fig. 6. This is a clear indication that the linear-response approximation is valid for this system in the classical limit, as was already discussed in previous work 12,53,55 . As a consequence, the classical limit of the aqueous ferrous-ferric system strongly resembles a spin-boson model, where the Wolynes theory is known to perform well. The same also explains the good agreement of Marcus theory with the exact classical rate, because Marcus theory, similarly to Wolynes theory, is known to perform well for this model. For a quantitative comparison, we note that the free-energy barrier according to Mar-  cus theory is Λ/4 = 27.1 mHartree, which is in excellent agreement with that found from Wolynes theory in the classical limit (27 ± 0.12 mHartree). Note that this is not always the case for asymmetric reactions, where the linear-response approximation is commonly seen to break down. 75 The aqueous ferrous-ferric system therefore does not exhibit a clear break-down of Wolynes theory in the classical limit, in contrast to the model systems tested in Ref. 40.
The investigation of the second qualitative indicator of breakdown of Wolynes theory is the distribution of values of σ = σ λ * (x) which are sampled in the unconstrained ensemble. If this distribution had a negligible population at σ = 0, which was the case for the system under study in Ref. 40, it would imply that the paths being sampled have no connection to the energy-conserving instanton paths and would be a clear sign of the break-down of Wolynes theory. However, as shown in Fig. 3, Wolynes calculations in the classical limit sample a uni-modal distribution peaked around σ = 0. This is similar to the picture that one would obtain for the spin-boson model for which Wolynes theory gives accurate rate predictions. This result therefore neither confirms nor disproves a break-down of Wolynes theory in this system.
A more detailed observation can be made from a comparison to another quantum rate theory, namely instanton theory. Instanton theory is not rigorously applicable to reactions in solution 25 and we cannot use it to calculate the rate. However, we can nonetheless acquire a qualitative insight into the tunnelling pathways of the system by obtaining a set of optimised instanton paths 27 on different solvent configurations randomly taken from a MD simulation (250 configurations). The water molecules beyond a radius of 5 Å of either Fe ion were fixed and only water molecules within this circumference were optimised in the instanton calculations (approximately 36 flexible water molecules). Analysing the ensemble of instantons clearly shows the multi-transition-state nature of the aqueous ferrous-ferric system. As shown in Fig. 7, the instantons have a range of different order parameters λ , which are found by a stationary-action principle. 23 A distribution of instantons occurs because even though the system is globally symmetric, it is locally asymmetric around each instanton and this is another sign that the system differs from the spin-boson model, which has only one instanton with λ = 0.5. The order parameters for each instanton cannot all be simultaneously satisfied by the single choice made by Wolynes theory of λ * = 0.5, which suggests that Wolynes theory may break down for this system and overpredict the rate. However, a broad uni-modal distribution as observed here is clearly a much safer scenario for Wolynes theory than the system tested in Ref. 40 for which this equivalent plot would have two peaks on either side of λ * . The fact that in this case the distribution is uni-modal and centred at λ = 0.5 (with a standard deviation of 0.07) suggests that the break down will be less severe.

Discussion of the GR-QTST result
GR-QTST, in contrast to Wolynes theory, can treat systems with multiple transition states correctly 40 , and similar to Wolynes theory, gives an accurate result for the spin-boson model also in the quantum limit 39 . In fact, the close agreement of GR-QTST and Wolynes theory for a spin-boson model 39 is another argument towards the aqueous ferrous-ferric system being badly approximated by a spin-boson model in the quantum limit, because for this system the two rate theories disagree (see Table 1).
Next to the discrepancy in rates, the curvature of the constrained free energy F c (λ ) shown in Fig. 8 can be utilised to argue against the applicability of linear-response theory and therefore the approximation of this system by a spin-boson model. Unlike the observation for a spin-boson model, for the atomistic system under study the constrained free energy F c (λ ) is curved upwards in extreme regions of the order parameter λ , i.e. far from the optimal order parameter λ * . In the earlier investigation of GR-QTST on spin-boson models 39 , we found that the constrained free energy F c (λ ) curves down when going away from the optimal order parameter λ * . This curvature behaviour becomes more prominent with an increased number of degrees of freedom and is already significant for 8 degrees of freedom, which can be understood based on the analysis presented in Ref. 39. The fact that the curvature behaviour differs from this is a further indicator for the aqueous ferrous-ferric system not being well described by a spin-boson model.
The curvature of the constrained free energy F c (λ ) is also of interest as it is an indication of the size-consistency error of GR-QTST. In Ref. 39 we showed that as more degrees of freedom are added to the system, the plot of F c (λ ) becomes more and more curved. However, we also argued that no matter how large the system, it becomes flat in the classical limit, and that for a spin-boson model, the value of F c (λ * ) remains stable even in the quantum case. There is however the possibility that this could lead to an error for more complex systems such as the aqueous ferrous-ferric reaction studied here. However, in the broad vicinity of the stationary point λ * the curvature of the constraint free energy is approximately flat. This is a good sign that there is no serious error being made by the GR-QTST method.

New and old controversies of the aqueous ferrous-ferric electron transfer
Already 30 years ago the investigation of the aqueous ferrousferric system led to a broad range of rate predictions, 33,53,55 and no conclusive argument could be made at the time for which was correct. We reopen this controversy by adding the results of our new GR-QTST approach to the discussion. When comparing the rate predictions to the earlier papers 33,53 , the enhancement due to tunnelling effects is of particular interest. This is quantified by the tunnelling enhancement factor Γ as introduced in Eqn. (16) and listed for the different rate theories employed in Table 2.
As reported in Table 2, our Wolynes rate calculations are in good agreement with the similar calculations performed by Bader et al. 33 . There are a couple of minor differences between the calculations besides the use of different water models. The rate Table 2 Tunnelling enhancement factors as defined in Eqn. (16) or accordingly. Note that the tunnelling enhancement factor reported by Bader et al. 33  predictions reported by Bader et al. only include the exponential contribution to the tunnelling enhancement factor as they assumed that the prefactor ratio would not be significant. From our Wolynes calculations, we find that the exponential contribution to the tunnelling factor is 58 ± 10 and the prefactor contribution is 1.2. Nonetheless, we find that these subtle differences do not distort the result significantly. The enhancement factor obtained from Wolynes theory however, is significantly higher than that predicted by Chandler and co-workers when employing the spin-boson model, which confirms our argument that the aqueous ferrous-ferric system cannot be well represented by a spin-boson model. In addition, the tunnelling enhancement factors from Wolynes theory and GR-QTST differ significantly, which is yet another argument towards the inapplicability of the spin-boson model here. In the quantum limit, the linear-response approximation therefore breaks down. Because of the large differences between Wolynes theory and GR-QTST and the fact that they predict very different tunnelling enhancement factors, no final statement on the exact contribution of quantum tunnelling to the reaction rate can be made.
Ultimately, no exact quantum-mechanical rate for the aqueous ferrous-ferric system can be obtained, but the comparison to experimental results is another possible aspect to investigate. For the ferrous-ferric system experimental isotope effects are available and the presence of a kinetic isotope effect is proof that nuclear quantum effects play a role in this reaction as the classical rate does not depend on the masses of the atoms. The experimental estimate of the isotope effect is in the range of k H exp /k D exp = 1.7 − 2.0 76,77 , which compares well to the ratio k H Wolynes /k D Wolynes = 2.0 ± 0.4 that we find by performing Wolynes calculations for both isotopes ((F u − F 0 ) H = 23.1 ± 0.11 mHartree and (F u − F 0 ) D = 23.8 ± 0.12 mHartree). Our prediction of the isotope effect from the GR-QTST calculations is k H /k D = 1.7 ± 0.3 and thus also lies within the range of experimental findings. The GR-QTST prediction is slightly lower than that of Wolynes theory, because the predicted tunnelling enhancement factor is smaller (see Table 2). Bader et al. 33 found an isotope factor of k H /k D = 2.6 ± 0.5. Our Wolynes results are accordingly within the error bars of the results obtained in the early papers 33 .
Although our main focus is on the comparison of different quantum rate theories, in order to justify our comparison with the experimental isotope effect, we must also consider the accuracy of the atomistic model. Marcus and co-workers predict a significantly lower tunnelling enhancement factor 55,56 than Chandler and co-workers using a spin-boson model but with a a different spectral density. They proposed in Ref. 55 that the discrepancy in rate predictions might be due to the neglect of polarisation and flexibility of water 33 . Calculations by Blumberger and coworkers 70,71 on a similar system (Ru 2+ -Ru 3+ ) however do not offer arguments for Marcus' hypothesis. Our results also suggest that Chandler's finding are robust in terms of the choice of a rigid water model as using a better flexible water model does not significantly change the tunnelling enhancement. However we do note that the experimentally estimated reorganisation energy (2.1 eV) 70,78 is lower than the one we predicted (2.95 eV), and a lower barrier could result in reduced tunnelling effects.
This work therefore aims to reopen the discussion on the question of tunnelling enhancement in the atomistic model of the ferrous-ferric electron transfer as presented here. We have shown that a number of different methods obtain a number of different rate predictions and currently no decisive argument for either can be made. This is reminiscent of the controversy surrounding the quantum tunnelling effect in the Azzouz-Borgis model of proton transfer, for which various approximate quantum rate theories do not agree [79][80][81][82] . Any further study of these problems however provides valuable insights into the accuracy, applicability and pitfalls of the quantum methods applied.

Conclusions
In this study, we have presented the first application of GR-QTST to an atomistic system of electron transfer and thereby obtained estimates for the reaction rate, isotope effects and tunnelling enhancement factors from this new method. The aqueous ferrousferric electron-transfer reaction has been extensively studied in the past and yet no conclusive quantitative answer for the contribution of nuclear tunnelling to the reaction rate has been given. In fact the previously predicted tunnelling enhancement factors span an order of magnitude 33,53,55,56 and our new prediction is at the lower end of this range.
All the methods tested reproduce the correct rate in the classical limit. GR-QTST is guaranteed to do this for any system, whereas Wolynes theory and Marcus theory are correct only if the linear-response approximation is valid. This implies that the classical limit of this reaction can be adequately described by a spin-boson model. This observation may lead one to believe that the spin-boson model is a feasible approximation of the aqueous ferrous-ferric system also in the quantum limit. However, we could show the unsuitability of this assumption by comparing Wolynes theory 31 and our newly developed GR-QTST 39,40 . In an earlier study it was observed that GR-QTST and Wolynes theory predict the same rate for a spin-boson model in the quantum limit 39 . In contrast to this, the two theories predict tunnelling enhancement factors that are different by a factor of 6 for the atomistic system, therefore making it impossible to argue that a spin-boson model is a good approximation in this case.
The obvious next question aims at resolving the disagreement of Wolynes theory and GR-QTST and in order to do so the possible pitfalls of each method were investigated. We have explained that there is a risk that Wolynes theory may break down and overpredict the rate, especially in liquid systems like that under study. This may occur whenever a reaction contains multiple distinct transition states 40 . We optimised a set of instantons in the system and found that their order parameters, λ , were distributed around λ * = 0.5. This implies that there are indeed multiple transition states in this system, although they are more similar to each other than the extreme cases studied in Ref. 40. The error made by Wolynes theory is thus expected to be less severe, but may still exist to some extent. GR-QTST is also not without its flaws and we have shown in previous work that the theory may suffer from size inconsistency. In a model system the addition of many degrees of freedom leads to a strongly λ -dependent F c (λ ) curve, which may degrade the rate predictions. Nonetheless, for the high-dimensional aqueous ferrous-ferric system, we observe only a minor curving of the constrained free energy of GR-QTST far from the optimal order parameter λ * , which is not expected to significantly degrade the result.
Considering the discrepancy in the predicted tunnelling enhancement factors makes it however apparent that at least one theory must be inaccurate for this system. A hypothesis that Wolynes theory is overpredicting the rate due to the multiinstanton nature of the system cannot be excluded from the results of our calculations. However, none of the studies provide us with a conclusive case to prove this statement, nor to rule out the possibility of an error on the part of GR-QTST.
By revisiting the controversial question of the nuclear tunnelling effect in this system, we could however show that the dynamics in this system deviates from those of the spin-boson model and raises the question of the applicability of this model for simulating atomistic systems. We hope to further elucidate the question of appropriateness of these quantum rate theories, by applying GR-QTST to an asymmetric atomistic system, where Wolynes theory may show more distinctive break-down behaviour.

Conflicts of interest
There are no conflicts to declare.