Rate constant expressions for vibronically nonadiabatic proton transfer and proton-coupled electron transfer reactions are presented and analyzed. The regimes covered include electronically adiabatic and nonadiabatic reactions, as well as high-frequency and low-frequency proton donor-acceptor vibrational modes. These rate constants differ from previous rate constants derived with the cumulant expansion approach in that the logarithmic expansion of the vibronic coupling in terms of the proton donor-acceptor distance includes a quadratic as well as a linear term. The analysis illustrates that inclusion of this quadratic term in the framework of the cumulant expansion framework may significantly impact the rate constants at high temperatures for proton transfer interfaces with soft proton donor-acceptor modes that are associated with small force constants and weak hydrogen bonds. The effects of the quadratic term may also become significant in these regimes when using the vibronic coupling expansion in conjunction with a thermal averaging procedure for calculating the rate constant. In this case, however, the expansion of the coupling can be avoided entirely by calculating the couplings explicitly for the range of proton donor-acceptor distances sampled. The effects of the quadratic term for weak hydrogen-bonding systems are less significant for more physically realistic models that prevent the sampling of unphysical short proton donor-acceptor distances. Additionally, the rigorous relation between the cumulant expansion and thermal averaging approaches is clarified. In particular, the cumulant expansion rate constant includes effects from dynamical interference between the proton donor-acceptor and solvent motions and becomes equivalent to the thermally averaged rate constant when these dynamical effects are neglected. This analysis identifies the regimes in which each rate constant expression is valid and thus will be important for future applications to proton transfer and proton-coupled electron transfer in chemical and biological processes.

## I. INTRODUCTION

Vibronic couplings play a key role in determining the rate constants for proton transfer (PT) and proton-coupled electron transfer (PCET). In the Golden Rule formulation, the rate constant is proportional to the square of the vibronic coupling at the environmental configurations corresponding to the degeneracy of the reactant and product vibronic states. In this context, the vibronic coupling is the Hamiltonian matrix element between the reactant and product diabatic electron-proton vibronic states. Typically, these vibronic couplings are assumed to be only weakly dependent on the configuration of the environment (i.e., the Condon approximation). In the case of PT or PCET reactions, however, the proton donor-acceptor distance coordinate represents an exception. In typical two-dimensional (2D) potentials describing the PT interface in terms of the proton coordinate *q* and the proton donor-acceptor distance *R*, the coordinates *q* and *R* are strongly coupled, resulting in the strong dependence of the PT barrier height along the proton coordinate *q* on the proton donor-acceptor distance *R*. A major consequence of this effect is the strong dependence of the vibronic coupling on the proton donor-acceptor distance.

The importance of the dependence of the vibronic coupling on the proton donor-acceptor vibrational motion has been recognized in the literature. The effects of the proton donor-acceptor motion were explicitly included in theoretical treaments of PT^{1–5} and PCET^{6–8} reactions. The dependence of the nonadiabatic vibronic coupling *V* on the proton donor-acceptor distance *R* is often represented in exponential form by expanding the natural logarithm of the coupling around the equilibrium proton donor-acceptor distance *R*_{0}. When only the terms up to second order in *δR* = *R* − *R*_{0} are retained, this dependence is expressed as

where

and

Using this expansion, the vibronic coupling can be written in the following exponential form with *V*_{0} ≡ *V*(*R*_{0}):

Several previous treatments^{1,2,4,5,8,9} assumed that the terms linear in *δR* are dominant and that the quadratic term can be neglected. This approximation was justified by agreement with the full quantum mechanical treatment of both the *q* and *R* coordinates in various model PCET systems for the electronically nonadiabatic case, where the vibronic coupling is proportional to the *R*-dependent overlap integral between the reactant and product proton vibrational wavefunctions.^{8,10} In other cases,^{3,11,12} however, the quadratic term in Eq. (4) could potentially have a significant effect on the rate constant, especially in the case of the so-called deep tunneling or electronically adiabatic but vibronically nonadiabatic regime, where the vibronic coupling is represented by half the tunneling splitting in a double-well PT potential. These effects are also relevant to thermal averaging treatments of PT and PCET, in which the rate constant is averaged over the equilibrium distribution of proton donor-acceptor distances.^{6} In this treatment, the classical distribution function over the proton donor-acceptor distances may enable sampling of the relatively short distances at which the linear approximation for the logarithm of the vibronic coupling is insufficient, and the quadratic term may become necessary.

In this paper, we systematically analyze the effects of the quadratic term on the nonadiabatic rate constants for PT and PCET in various regimes. The aim of this study is to obtain deeper insight into how the characteristics of the PT interface, including the features of the potential energy surface and the effective mass and frequency of the proton donor-acceptor vibrational mode, impact the significance of the quadratic term. Section II summarizes the theoretical treatments for nonadiabatic PT and PCET reactions and provides expressions for the vibronic coupling in the relevant regimes. Section III presents results for a range of model systems, and Section IV summarizes the conclusions of this study.

## II. THEORY

The nonadiabatic theories of PT and PCET are based on the general Golden Rule formulation. In the conventional formulation, nonadiabatic transitions occur between the reactant and product quantum states interacting with the fluctuating environment in thermal equilibrium at a given temperature. The thermal fluctuations of the environment bring a pair of quantum states into degeneracy, thereby enabling a quantum transition with a probability that is proportional to the square of the quantum mechanical coupling between the corresponding states. Defining the reactant (I) and product (II) channels with Hamiltonians $H\u02c6I$ and $H\u02c6II$, which give rise to the manifolds of reactant and product vibronic states for the entire system, the Golden Rule rate constant expression can be written as a thermodynamic average over the reactant states (i.e., eigenstates of the Hamiltonian $H\u02c6I$) of the entire system:^{8}

The time-dependent quantities in the above expression are op erators for the coupling $V\u02c6$ and the energy gap $\Delta H\u02c6\u2261H\u02c6II\u2212H\u02c6I$ in the Heisenberg representation:

_{(−)}denotes a negatively time-ordered exponential, and $\cdots {I}\u2261Trexp\u2212\beta H\u02c6I\cdots /Trexp\u2212\beta H\u02c6I$, where

*β*

^{−1}=

*k*

_{B}

*T*, denotes a quantum mechanical thermal average over initial states.

The reactant and product vibronic states form two sets of quantum states in the electronically diabatic representation. As in Marcus-Levich-Dogonadze nonadiabatic ET theory,^{13,14} the diabatic electronic states are characterized by the reactant and product electronic charge distributions in the reaction complex. These diabatic electronic states are defined either by the transferring electron localized on the electron donor or acceptor species, as in the case of ET or PCET, or by the specific bonding pattern describing the proton bonded to the proton donor or acceptor, as in the case of PT. For PCET and PT, the proton donor-acceptor distance *R* warrants special treatment due to its special role in modulating the couplings between the reactant and product vibronic states. In the remainder of this section, we consider different treatments of this coordinate, emphasizing the effects of the linear and quadratic terms in the logarithmic expansion of the coupling matrix element given in Eq. (4). In the general rate constant expressions, we do not specify the form of the vibronic coupling, assuming that any coupling matrix element can be represented in the form of Eq. (4) near the equilibrium value of the proton donor-acceptor distance. Specific expressions for the vibronic coupling in the electronically adiabatic and nonadiabatic limits, as well as the intermediate regime, are provided at the end of Section II A.

### A. Reactant and product vibronic states and vibronic couplings

The reactant and product vibronic states are the eigenstates of the reactant and product channel Hamiltonians $H\u02c6I$ and $H\u02c6II$, respectively. Within each of these channels, the reactant and product states can be represented in the following form:

Here, $\Phi \mu J$ are the diabatic electron-proton vibronic wavefunctions parametrically dependent on the donor-acceptor distance *R* and the environmental coordinates *ξ* for the reactant (*J* = I) and product (*J* = II) vibronic manifolds with corresponding diabatic electronic charge distributions. The corresponding vibronic energies, $E\mu J(R,\xi )$, serve as vibronic potential energy surfaces for the remaining degrees of freedom, and the final quantization along *R* and *ξ* provides the vibrational wavefunctions of the environment, $\chi \mu mJ(R,\xi )$. In the linear response approximation, the vibronic surfaces $E\mu J(R,\xi )$ are represented as harmonic surfaces in the multidimensional space spanned by the proton donor-acceptor coordinate *R* and the environmental modes.

The matrix elements of the coupling operator $V\u02c6$ between the reactant and product states defined in Eq. (7) are given by

where

is the electron-proton vibronic coupling appearing in Eq. (4). According to the Condon approximation, this coupling is assumed to be independent of the environmental coordinates, but its dependence on the proton donor-acceptor distance cannot be neglected.

In the electronically diabatic channel representation, the hierarchy of adiabatic approximations between the electrons and the proton, as well as between the proton and the slower nuclear degrees of freedom within each channel, is usually well-justified. The adiabatic separation between the proton coordinate *q* and the proton donor-acceptor coordinate *R* is less obvious because these coordinates are strongly coupled in typical two-dimensional potentials describing hydrogen-bonded PT interfaces. However, previous work^{15} has established that this adiabatic separation, in which the proton coordinate *q* is treated as a “fast” coordinate and the proton donor-acceptor coordinate *R* is treated as a “slow” coordinate, is well-justified not only for electronically diabatic but also for a wide range of electronically adiabatic potentials, including double well potentials with high barriers (i.e., in the deep tunneling regime).

The vibronic coupling *V*(*R*) can be calculated from the proton potential energy curves associated with the reactant and product diabatic electronic states, as depicted in Figure 1(a), in conjunction with the electronic coupling *V*^{el} between these two states. Georgievskii and Stuchebrukhov^{16} derived a semiclassical expression for the general vibronic coupling *V*^{(sc)},

where *V*^{(ad)} = Δ^{(tun)}/2 is the vibronic coupling in the electronically adiabatic limit, equal to the half the tunneling splitting between the proton vibrational levels in the ground state adiabatic electronic potential (Figure 1(b)), and

In this expression, Γ(*x*) is the gamma function, and *p* = *τ*_{p}/*τ*_{e} is the electron-proton adiabaticity parameter defined in terms of the effective electronic transition time,

and the effective proton tunneling time,

Here, *V*^{el} is the electronic coupling between the reactant and product diabatic electronic states,^{17} |Δ*F*| is the difference between the slopes of the electronically diabatic potentials at the crossing point, and the tunneling velocity is

where *V*^{†} is the energy at which the potential energy curves cross, *E* is the tunneling energy associated with the unperturbed degenerate proton vibrational levels in the reactant and product diabatic potentials, and *m* is the mass of the proton.

The semiclassical expression for the vibronic coupling given in Eq. (10) spans the electronically adiabatic and nonadiabatic regimes. In the electronically adiabatic limit, *p* ≫ 1, *κ* = 1, and the vibronic coupling simplifies to *V*^{(ad)}. Several semiclassical expressions in terms of parameters associated with proton double well potentials have been derived for *V*^{(ad)}.^{5,11,12,18,19} In the electronically nonadiabatic limit, *p* ≪ 1, $\kappa =2\pi p1/2$, and the vibronic coupling reduces to

where *φ*^{(I)} and *φ*^{(II)} are the proton vibrational wavefunctions in the reactant and product wells, respectively (Figure 1(a)).

From a physical perspective, the electronically adia batic limit corresponds to the case in which the electrons respond instantaneously to the proton motion, resulting in the reaction proceeding entirely on the ground adiabatic electronic state. In contrast, the electronically nonadiabatic limit corresponds to the case in which the electrons are unable to respond instantaneously to the proton motion, thereby leading to involvement of the excited electronic state. Most PT and hydrogen atom transfer (HAT) reactions are electronically adiabatic, occurring on the ground adiabatic electronic state. In contrast, PCET reactions that involve concerted electron and proton transfer between different donors and acceptors (i.e., different orbitals, atoms, or groups) are typically electronically nonadiabatic. The distinctions between these types of reactions have been discussed extensively elsewhere.^{20–24}

### B. Cumulant expansion approach to Golden Rule nonadiabatic rate constant

In the approaches used previously for PT^{2,5} and PCET^{8} reactions in solution, the proton donor-acceptor coordinate *R* is assumed to be adiabatically separated from the electron and proton coordinates. Here, we limit our discussion to a single pair of reactant and product electron-proton vibronic states $\Phi \mu I$ and $\Phi \nu II$ to avoid complicated notation; the total rate constant is obtained by averaging over the reactant states and summing over the product states.

Expanding the energy gap in a Taylor series in the coordinate *R* up to the linear term and substituting the expansion in Eq. (4) for the vibronic coupling, we obtain the following expression for the reaction flux correlation function, *J*(*t*), in Eq. (5):^{5}

where $\delta E\u02c6(\tau )$ and $\delta R\u02c6(\tau )$ are the Heisenberg operators describing the time evolution of fluctuations of the energy gap $E\u02c6\u2261\Delta H\u02c6$ and proton donor-acceptor distance *R*, respectively. Moreover, $\alpha \u0303=\u2212i\u0127\alpha $, $\gamma \u0303=\u2212i\u0127\gamma $, and $\Lambda \u0303R$ is the derivative of the energy gap with respect to the proton donor-acceptor coordinate *R* evaluated at the minimum of the reactant vibronic energy surface. In the harmonic approximation, $\Lambda \u0303R$ is proportional to Δ*R*, the shift of the equilibrium value of the proton donor-acceptor distance in the product state relative to that in the reactant state. Note that for symmetric systems, the reactant and product equilibrium proton donor-acceptor distances are equal, and thus, $\Lambda \u0303R=0$.

Expanding the exponentials in Eq. (16) and collecting terms of the same order in *i*/ℏ, we obtain a series that can be represented as a cumulant expansion of the time-ordered exponential for the time-dependent variable $X\u02c6(\tau )$,

where

and double angular brackets denote the cumulants of $X\u02c6(\tau )$, which can be expressed in terms of the central moments and correlation functions of the time-dependent fluctuations of the energy gap and proton donor-acceptor distance, $\delta E\u02c6$ and $\delta R\u02c6$, respectively. These moments can in turn be expressed in terms of the corresponding cumulants. Next, we assume that these time-dependent fluctuations can be described by Gaussian processes corresponding to the time evolution on harmonic surfaces. For the proton donor-acceptor coordinate *R*, the Gaussian process $\delta R\u02c6(t)$ corresponds to the time evolution of a harmonic oscillator with an effective mass *M* and frequency Ω. In practice, the effective mass and frequency can be defined in terms of the normal modes of the system that have nonvanishing projections on the proton donor-acceptor axis and thus contribute to the fluctuations of the coordinate *R*.^{25} In addition, we assume that the corresponding harmonic potential along *R* is uncoupled from the other environmental degrees of freedom, which can be achieved by considering the normal modes of the entire system in the reactant state.

For Gaussian processes $\delta E\u02c6(\tau )$ and $\delta R\u02c6(\tau )$, all of the corresponding cumulants higher than second order vanish exactly. However, the cumulant expansion for the rate constant based on the expansion of the time-ordered exponential with the variable $X\u02c6(\tau )$ does not truncate even for Gaussian $\delta E\u02c6(\tau )$ and $\delta R\u02c6(\tau )$ due to the presence of the quadratic terms $\delta R\u02c62(t)$ in Eq. (18). The final expression for the reaction flux can be written as

where

and the operator notation is omitted for simplicity. Here, Δ*G* and Λ are the reaction free energy and the total reorganization energy, respectively, and the correlation functions are defined as

*n*th order polynomials in

*C*(

_{R}*t*) and $\delta R2=CR(0)$. This additional term arises because of the quadratic term in the logarithmic expansion of the coupling. The explicit expressions for $Kn(1)(t)$ and $Kn(2)(t)$ for

*n*= 1–8 (i.e., up to eighth order in

*γ*) are given in Table I.

n
. | $Kn(1)(t)$ . | $Kn+1(2)(t)$ . |
---|---|---|

1 | κ_{2,0} | $\kappa 1,12+2\kappa 2,0\kappa 1,1+\kappa 2,02$ |

2 | $\kappa 1,12+\kappa 2,02$ | $\kappa 1,13+3\kappa 2,0\kappa 1,12+3\kappa 2,02\kappa 1,1+\kappa 2,03$ |

3 | $\kappa 2,03+3\kappa 1,12\kappa 2,0$ | $\kappa 1,14+4\kappa 2,0\kappa 1,13+6\kappa 2,02\kappa 1,12+4\kappa 2,03\kappa 1,1+\kappa 2,04$ |

4 | $\kappa 1,14+6\kappa 2,02\kappa 1,12+\kappa 2,04$ | $\kappa 1,15+5\kappa 2,0\kappa 1,14+10\kappa 2,02\kappa 1,13+10\kappa 2,03\kappa 1,12+5\kappa 2,04\kappa 1,1+\kappa 2,05$ |

5 | $\kappa 2,05+10\kappa 1,12\kappa 2,03+5\kappa 1,14\kappa 2,0$ | $\kappa 1,16+6\kappa 2,0\kappa 1,15+15\kappa 2,02\kappa 1,14+20\kappa 2,03\kappa 1,13+15\kappa 2,04\kappa 1,12+6\kappa 2,05\kappa 1,1+\kappa 2,06$ |

6 | $\kappa 1,16+15\kappa 2,02\kappa 1,14+15\kappa 2,04\kappa 1,12+\kappa 2,06$ | $\kappa 1,17+7\kappa 2,0\kappa 1,16+21\kappa 2,02\kappa 1,15+35\kappa 2,03\kappa 1,14+35\kappa 2,04\kappa 1,13+21\kappa 2,05\kappa 1,12+7\kappa 2,06\kappa 1,1+\kappa 2,07$ |

7 | $\kappa 2,07+21\kappa 1,12\kappa 2,05+35\kappa 1,14\kappa 2,03+7\kappa 1,16\kappa 2,0$ | $\kappa 1,18+8\kappa 2,0\kappa 1,17+28\kappa 2,02\kappa 1,16+56\kappa 2,03\kappa 1,15+70\kappa 2,04\kappa 1,14+56\kappa 2,05\kappa 1,13+28\kappa 2,06\kappa 1,12+8\kappa 2,07\kappa 1,1+\kappa 2,08$ |

8 | $\kappa 1,18+28\kappa 2,02\kappa 1,16+70\kappa 2,04\kappa 1,14+28\kappa 2,06\kappa 1,12+\kappa 2,08$ | $\kappa 1,19+9\kappa 2,0\kappa 1,18+36\kappa 2,02\kappa 1,17+84\kappa 2,03\kappa 1,16+126\kappa 2,04\kappa 1,15+126\kappa 2,05\kappa 1,14+84\kappa 2,06\kappa 1,13+36\kappa 2,07\kappa 1,12+9\kappa 2,08\kappa 1,1+\kappa 2,09$ |

n
. | $Kn(1)(t)$ . | $Kn+1(2)(t)$ . |
---|---|---|

1 | κ_{2,0} | $\kappa 1,12+2\kappa 2,0\kappa 1,1+\kappa 2,02$ |

2 | $\kappa 1,12+\kappa 2,02$ | $\kappa 1,13+3\kappa 2,0\kappa 1,12+3\kappa 2,02\kappa 1,1+\kappa 2,03$ |

3 | $\kappa 2,03+3\kappa 1,12\kappa 2,0$ | $\kappa 1,14+4\kappa 2,0\kappa 1,13+6\kappa 2,02\kappa 1,12+4\kappa 2,03\kappa 1,1+\kappa 2,04$ |

4 | $\kappa 1,14+6\kappa 2,02\kappa 1,12+\kappa 2,04$ | $\kappa 1,15+5\kappa 2,0\kappa 1,14+10\kappa 2,02\kappa 1,13+10\kappa 2,03\kappa 1,12+5\kappa 2,04\kappa 1,1+\kappa 2,05$ |

5 | $\kappa 2,05+10\kappa 1,12\kappa 2,03+5\kappa 1,14\kappa 2,0$ | $\kappa 1,16+6\kappa 2,0\kappa 1,15+15\kappa 2,02\kappa 1,14+20\kappa 2,03\kappa 1,13+15\kappa 2,04\kappa 1,12+6\kappa 2,05\kappa 1,1+\kappa 2,06$ |

6 | $\kappa 1,16+15\kappa 2,02\kappa 1,14+15\kappa 2,04\kappa 1,12+\kappa 2,06$ | $\kappa 1,17+7\kappa 2,0\kappa 1,16+21\kappa 2,02\kappa 1,15+35\kappa 2,03\kappa 1,14+35\kappa 2,04\kappa 1,13+21\kappa 2,05\kappa 1,12+7\kappa 2,06\kappa 1,1+\kappa 2,07$ |

7 | $\kappa 2,07+21\kappa 1,12\kappa 2,05+35\kappa 1,14\kappa 2,03+7\kappa 1,16\kappa 2,0$ | $\kappa 1,18+8\kappa 2,0\kappa 1,17+28\kappa 2,02\kappa 1,16+56\kappa 2,03\kappa 1,15+70\kappa 2,04\kappa 1,14+56\kappa 2,05\kappa 1,13+28\kappa 2,06\kappa 1,12+8\kappa 2,07\kappa 1,1+\kappa 2,08$ |

8 | $\kappa 1,18+28\kappa 2,02\kappa 1,16+70\kappa 2,04\kappa 1,14+28\kappa 2,06\kappa 1,12+\kappa 2,08$ | $\kappa 1,19+9\kappa 2,0\kappa 1,18+36\kappa 2,02\kappa 1,17+84\kappa 2,03\kappa 1,16+126\kappa 2,04\kappa 1,15+126\kappa 2,05\kappa 1,14+84\kappa 2,06\kappa 1,13+36\kappa 2,07\kappa 1,12+9\kappa 2,08\kappa 1,1+\kappa 2,09$ |

In the high-temperature approximation for the solvent fluctuations and for a simple Debye model of solvent relaxation, the correlation function of the energy gap becomes a single exponential.^{26} In this case, the double time integral with the energy gap correlation function in Eq. (19) can be approximated by a Gaussian,

where *λ*_{s} is the solvent reorganization energy. Note that the simple Debye model of solvent relaxation breaks down when the short time scale solvent dynamics is important, in which case the energy gap correlation function may include Gaussian components to describe the short time scale relaxation processes.^{27,28} Using the readily available analytical expressions for the variance 〈*δR*^{2}〉 and correlation function $\delta R(0)\u2009\delta R(\tau )$ of the quantum harmonic oscillator with mass *M* and frequency Ω,^{29} we obtain the following general expression for the reaction flux *J*(*t*):

In this expression, $\lambda \alpha =\u01272\alpha 22M$ and $\lambda \gamma =\u01272\gamma 2M$ are the coupling reorganization energies corresponding to the linear and quadratic terms, respectively, in the logarithmic expansion of the coupling, $\zeta =coth\beta \u0127\Omega /2$, and $K\u0303n(i)=Kn(i)/\delta R2n$ for *i* = 1, 2. The remaining parameters are defined exactly as in our previous work:^{8}

The additional exponential factor arising from the quadratic term in the coupling expansion is an infinite sign-alternating series that in certain regimes may cause convergence difficulties for the calculation of the reaction flux given in Eq. (23). The key temperature-dependent convergence parameter is *η*_{γ} = *λ*_{γ}*ζ*/*ħ*Ω, which depends on *γ* as well as on the mass *M* and the frequency Ω of the proton donor-acceptor mode. In the low-frequency limit for the proton donor-acceptor mode, *βħ*Ω ≪ 1, the convergence parameter becomes $\eta \gamma =\gamma /\beta M\Omega 2$. Thus, at high temperatures and for a very soft proton donor-acceptor mode with small force constant *M*Ω^{2}, the convergence may become slow. The convergence should improve significantly in the high-frequency limit, *βħ*Ω ≫ 1, where $\eta \gamma =\lambda \gamma /\u0127\Omega \u226a1$.

For the case in which only the linear term in the logarithmic coupling expansion is retained (i.e., *γ* = 0), closed analytical forms of the rate constant have been derived in the low- and high-frequency regimes for the proton donor-acceptor mode.^{5,8} In the case when the quadratic term is included (i.e., *γ*≠0), a closed analytical expression for the rate constant can be obtained if we neglect the equilibrium dynamics of the proton donor-acceptor coordinate *R* by replacing the time correlation function $CR(t)=\delta R(0)\u2009\delta R(t)$ by its value at *t* = 0, $CR(0)=\delta R2$. This approximation for the proton donor-acceptor mode is valid when the characteristic time scale of the correlation function *C _{R}*(

*t*) is much longer than the time scale associated with the solvent damping term given in Eq. (22). In previous work,

^{7}we analyzed the time dependence of the correlation function

*C*(

_{R}*t*) calculated from molecular dynamics simulations of model PCET systems and its effects on the rate constant for the linear logarithmic coupling expansion. For typical PT interfaces, we found that the approximation

*C*(

_{R}*t*) ≈

*C*(0) is valid at ambient temperatures for solvents with large reorganization energy. In this case, the infinite series in Eq. (23) reduces to a time-independent prefactor, and the time integral of the reaction flux correlation function can be readily evaluated. For PT interfaces with $\Lambda \u0303R=0$, the rate constant acquires the following analytical form:

_{R} The low-frequency (*βħ* Ω ≪ 1) and high-frequency (*βħ*Ω ≫ 1) limits for the proton donor-acceptor mode can be obtained from the above expression by substituting the corresponding limits for $\zeta =coth\beta \u0127\Omega /2$.

Unfortunately, in the general case, when the time dependence of *C _{R}*(

*t*) is retained and the interference effects between the short time dynamics of the proton donor-acceptor mode and the fluctuations of the energy gap are included, the inclusion of the quadratic term in the coupling expansion leads to expressions that are not easily simplified into this type of analytical form. In this case, the effects of the quadratic term on the rate constant cannot be represented as a simple prefactor in the rate constant expression. In our calculations of the nonadiabatic rate constants for model systems discussed below, we use Eq. (5) with the full reaction flux given in Eq. (23) and perform the integration numerically.

### C. Thermal averaging over proton donor-acceptor motion

The alternative thermal averaging approach^{6} assumes that the proton donor-acceptor mode remains in thermal equilibrium during the reaction and that the nonadiabatic transitions or proton tunneling events occur independently at all values of *R*. This assumption is valid when the relaxation along *R* is fast compared to the rate of nonadiabatic transitions and is similar to the slow reaction limit in the Sumi-Marcus model for electron transfer.^{30} In this case, ignoring the weak dependence of the reorganization energy and the reaction free energy on *R*, the nonadiabatic rate constant is given by the following expres sion:

where the first prefactor is obtained by simple thermal averaging of the squared vibronic coupling over an equilibrium distribution of the proton donor-acceptor distances. For the exponential form of the coupling given in Eq. (4), the average is expressed as

where *P*(*R*) is the equilibrium distribution function for proton donor-acceptor distances.

In the low-frequency (*βħ*Ω ≪ 1) and high-frequency (*βħ*Ω ≫ 1) limits for a harmonic oscillator representing the proton donor-acceptor vibrational mode, the integral in Eq. (27) can be evaluated analytically, resulting in the following analytical expressions:

As pointed out previously,^{31} in the absence of the quadratic term in the coupling expansion, the averaging of the squared coupling with the classical harmonic distribution function leads to a prefactor that coincides with the prefactor in the rate constant expression derived with the cumulant expansion approach in the low-frequency limit for the proton donor-acceptor mode. When the quadratic term is included in the coupling expansion, the expressions given in Eqs. (28) and (29) coincide with the low-frequency (*βħ*Ω ≪ 1) and high-frequency (*βħ*Ω ≫ 1) limits of the prefactor in Eq. (25). Thus, the thermal averaging of the rate constant with the equilibrium distribution function for the proton donor-acceptor mode, resulting in the rate constant expression in Eq. (26), corresponds to an approximation that completely neglects the effects from the dynamical interference between the proton donor-acceptor motion and the fluctuations of the energy gap. This relation between the cumulant expansion and thermally averaged rate constant expressions can be proven for any form of the vibronic coupling. From a theoretical perspective, the cumulant expansion approach is more general and consistent because the equilibrium dynamics along the *R* coordinate is treated explicitly on the same level as the dynamics of the energy gap.

The expression in Eq. (28) is similar to the result obtained by Trakhtenberg *et al.*^{1} for hydrogen tunneling reactions with a single promoting mode (cf. Eq. (8) in Ref. 1). Analyzing the tunneling splittings in proton potentials derived from Morse functions for chemical bonds, the authors concluded that the additional prefactors arising from the quadratic term in the logarithmic coupling expansion can be safely neglected. However, other studies^{3,12,32} have found that at high temperatures and for soft PT interfaces with very small force constants *M*Ω^{2}, the second exponential factor in Eq. (28) can become small and thus can significantly reduce the rate constant. When only the linear term in the coupling expansion is included, the coupling is overestimated at short proton donor-acceptor distances, and for soft PT interfaces at high temperatures, the distribution function *P*(*R*) is broad enough to allow sampling of these short distances. On the other hand, the use of more realistic potentials along the coordinate *R* for hydrogen-bonded PT interfaces, such as Morse-like potentials with repulsive walls preventing the proton donor and acceptor from unphysical close contact, reduces the effects of the quadratic term due to decreased sampling probability at unphysical short donor-acceptor distances.

In an alternative, explicit thermal averaging approach, the vibronic coupling is calculated analytically or numerically at each proton donor-acceptor distance along a grid, followed by numerical evaluation of the first integral in Eq. (27).^{25,33–36} Thus, the vibronic couplings are calculated explicitly for the entire range of proton donor-acceptor distances sampled. This approach does not rely on the expansion of the vibronic coupling near the equilibrium proton donor-acceptor distance and therefore avoids any assumption of linear or quadratic behavior.

### D. Two-dimensional quantum treatment of proton and donor-acceptor coordinates

Another approach, which should be valid for a wide range of proton donor-acceptor mode effective masses *M* and frequencies Ω, employs a 2D quantum treatment of the PT interface. In this case, the *R* coordinate is treated on the same level as the electron and proton coordinates, and the reactant and product vibronic states are defined as solutions of the vibrational Schrödinger equation for the proton and its donor and acceptor moving on 2D diabatic electronic potentials *U ^{J}*(

*q*,

*R*) corresponding to the reactant and product diabatic electronic states,

*J*= I, II. The shape of these diabatic potentials is assumed to be independent of the environmental coordinates

*ξ*according to the definition of diabatic electronic states. In the case of 2D harmonic potentials

*U*(

^{J}*q*,

*R*), the solutions of the vibrational Schrödinger equation can be obtained analytically,

^{12}and the vibronic couplings can be evaluated analytically in the electronically nonadiabatic regime when they are proportional to the overlap integrals between the 2D vibrational wavefunctions, $\psi \mu J(q,R)$, of the reactant and product states. A detailed analysis of such overlap integrals has been performed recently.

^{12}For more realistic anharmonic 2D potentials, the 2D vibrational wavefunctions $\psi \mu J(q,R)$ can be calculated numerically using, for example, Fourier Grid Hamiltonian (FGH) methods.

^{37,38}

Alternatively, the adiabatic separation between the proton coordinate *q* and the donor-acceptor coordinate *R* can be invoked by treating *q* as a fast coordinate and *R* as a slow coordinate. In this case, the reactant and product vibronic states are products of the adiabatic proton vibrational wavefunctions $\chi \mu J(q|R)$, parametrically dependent on the *R* coordinate, and vibrational wavefunctions for the *R* coordinate $\varphi \mu mJ(R)$. As mentioned in Section II A, such an approximation is justified for hydrogen-bonded PT interfaces due to the specific form of the coupling between the coordinates *q* and *R* and the large difference in associated masses.^{15} For adiabatically separated *q* and *R* in the reactant and product potentials, the nonadiabatic rate constant for a Debye solvent is given by

where $P\mu mI$ is the Boltzmann probability for the reactant state $\chi \mu I\varphi \mu mI$ and Δ*G*_{μm,νn} is the reaction free energy for the pair of vibronic states *μm* and *νn*.

In previous work,^{8} we compared the full 2D quantum treatment of coordinates *q* and *R* with other approaches for a series of models with anharmonic diabatic electronic potential energy surfaces. Numerical calculations of the nonadiabatic rate constants in the electronically nonadiabatic regime for various models illustrated that at ambient temperatures of ∼300 K and typical *R*-mode frequencies of 100–500 cm^{−1} and masses of 10–100 amu, the full 2D quantum and explicit thermal averaging approaches provide very similar results. Note that the general cumulant expansion approach described in Section II B ^{2,5,8,10} also treats both the *q* and *R* coordinates quantum mechanically when the quantum mechanical correlation function for the harmonic *R*-mode is utilized. Moreover, similar to the thermal averaging approach, the 2D quantum approach neglects the interference effects arising from the dynamics of the proton donor-acceptor motion that are included in the cumulant expansion approach.

## III. MODEL CALCULATIONS

In this section, we analyze the effects of the quadratic term in the logarithmic expansion of the vibronic coupling on the rate constants for model systems spanning various regimes.

### A. Model potentials

For the analysis in the electronically nonadiabatic regime, we employed two simple analytical models with harmonic and Morse proton potentials, denoted 1D-HARM and 1D-MORSE, respectively. These simple symmetric models assume that the shape of the diabatic proton potentials does not depend on the proton donor-acceptor coordinate *R*, which affects only the distance between the minima on the reactant and product potentials. The motion along *R* is assumed to be harmonic with effective mass *M* and frequency Ω. The harmonic proton potentials in the 1D-HARM model have an equilibrium distance of 1.0 Å, proton mass of 1.007 amu, and frequency of 2900 cm^{−1}. The parameters of the Morse proton potentials in the 1D-MORSE model are similar to the parameters used to describe the nonadiabatic PCET reaction in soybean lypoxygenase:^{39} the bond dissociation energy *D* = 80 kcal/mol, the Morse parameter *β* = 2.1189 Å^{−1}, and the bond equilibrium distance *r*^{0} = 1.0 Å. The proton vibrational wavefunctions and overlap integrals for the harmonic and Morse proton potentials were calculated analytically,^{40,41} and the coupling attenuation parameters *α* and *γ* were calculated by evaluating the derivatives in Eqs. (2) and (3) numerically.

For the analysis in the intermediate and electronically adiabatic regimes, we used 2D hybrid non-harmonic proton potentials that combine a London-Eyring-Polanyi-Sato (LEPS) potential *U*^{LEPS}(*q*, *R*) for triatomic collinear systems and a harmonic potential *U*^{harm}(*R*) along the proton donor-acceptor coordinate *R*. The electronically diabatic potentials are represented by the diagonal elements of the 2 × 2 matrix of the electronic Hamiltonian with matrix elements:

where *V*^{el} is a constant electronic coupling parameter chosen to be 10 kcal/mol. The corresponding adiabatic potential $U0(ad)(q,R)$ for the ground electronic state is obtained by diagonalization of this [2 × 2] matrix.

This potential was used in one of our earlier model studies,^{8} and the reader is referred to Ref. 8 for technical details and the values of the parameters for the LEPS potentials. The two symmetric models are denoted as LEPS-M10-Ω200 and LEPS-M50-Ω200 according to the values of the effective mass *M* and frequency Ω of the donor-acceptor mode. The contour plots for the diabatic 2D reactant and product potentials are shown in Figure 2 for one of the models. We point out that these 2D model potentials include the effects of Duschinski rotations that have been considered to be important.^{12} We used the FGH method for the numerical calculation of the tunneling splittings in the electronically adiabatic proton potential for one-dimensional slices along the proton coordinate *q*. From these tunneling splittings, we obtained the vibronic coupling in the electronically adiabatic limit and calculated the semiclassical vibronic couplings in the intermediate regime using Eq. (10). The corresponding coupling attenuation parameters *α* and *γ* were calculated numerically using Eqs. (2) and (3), respectively.

### B. Analysis in electronically nonadiabatic regime

In the electronically nonadiabatic regime, the vibronic coupling is of the form given in Eq. (15), which is proportional to the overlap integral between the reactant and product proton vibrational wavefunctions associated with the reactant and product diabatic electronic states, respectively. To estimate the significance of the quadratic term in the logarithmic coupling expansion, we calculated the linear and quadratic attenuation parameters, *α* and *γ*, respectively, for the proton vibrational overlap integrals for the model systems with harmonic and Morse proton potentials. The calculated values of these attenuation parameters for proton donor-acceptor distances ranging from 2.7 to 3.0 Å are provided in Table II. Because of the Gaussian nature of the harmonic overlap integrals, the *γ* attenuation parameters are larger for the harmonic potentials than for the Morse potentials. The asymmetry of the Morse potentials leads to slightly enhanced proton delocalization toward the acceptor site, leading to smaller values of the *γ* parameters.

R_{0}, Å
. | α (H)
. | γ (H)
. | α (D)
. | γ (D)
. | λ_{α} (H)
. | λ_{γ} (H)
. | λ_{α} (D)
. | λ_{γ} (D)
. |
---|---|---|---|---|---|---|---|---|

1D-HARM model | ||||||||

2.7 | 30.32 | 43.32 | 42.87 | 61.25 | 4.432 | 0.209 | 8.860 | 0.295 |

2.8 | 34.67 | 43.32 | 48.99 | 61.25 | 5.789 | 0.209 | 11.57 | 0.295 |

2.9 | 38.99 | 43.32 | 55.12 | 61.25 | 7.327 | 0.209 | 14.65 | 0.295 |

3.0 | 43.32 | 43.32 | 61.24 | 61.25 | 9.045 | 0.209 | 18.08 | 0.295 |

1D-MORSE model | ||||||||

2.7 | 19.83 | 20.64 | 28.69 | 29.18 | 0.379 | 0.0199 | 0.79 | 0.0281 |

2.8 | 21.79 | 18.57 | 31.46 | 26.25 | 0.458 | 0.0179 | 0.95 | 0.0253 |

2.9 | 23.55 | 16.70 | 33.95 | 23.61 | 0.535 | 0.0161 | 1.11 | 0.0228 |

3.0 | 25.14 | 15.03 | 36.19 | 21.24 | 0.609 | 0.0145 | 1.26 | 0.0205 |

R_{0}, Å
. | α (H)
. | γ (H)
. | α (D)
. | γ (D)
. | λ_{α} (H)
. | λ_{γ} (H)
. | λ_{α} (D)
. | λ_{γ} (D)
. |
---|---|---|---|---|---|---|---|---|

1D-HARM model | ||||||||

2.7 | 30.32 | 43.32 | 42.87 | 61.25 | 4.432 | 0.209 | 8.860 | 0.295 |

2.8 | 34.67 | 43.32 | 48.99 | 61.25 | 5.789 | 0.209 | 11.57 | 0.295 |

2.9 | 38.99 | 43.32 | 55.12 | 61.25 | 7.327 | 0.209 | 14.65 | 0.295 |

3.0 | 43.32 | 43.32 | 61.24 | 61.25 | 9.045 | 0.209 | 18.08 | 0.295 |

1D-MORSE model | ||||||||

2.7 | 19.83 | 20.64 | 28.69 | 29.18 | 0.379 | 0.0199 | 0.79 | 0.0281 |

2.8 | 21.79 | 18.57 | 31.46 | 26.25 | 0.458 | 0.0179 | 0.95 | 0.0253 |

2.9 | 23.55 | 16.70 | 33.95 | 23.61 | 0.535 | 0.0161 | 1.11 | 0.0228 |

3.0 | 25.14 | 15.03 | 36.19 | 21.24 | 0.609 | 0.0145 | 1.26 | 0.0205 |

The important quantities affecting the nonadiabatic rate constant obtained by the cumulant expansion approach, based on the reaction flux given in Eq. (23), are the coupling reorganization energies *λ*_{α} and *λ*_{γ}, which are provided in Table II. The values of *λ*_{γ} are an order of magnitude smaller than the already very small values of *λ*_{α} for *M* = 10 amu, and they become even smaller as the effective mass of the proton donor-acceptor mode increases. The small magnitude of *λ*_{γ} results in small values of the convergence parameter *η*_{γ} that ensures convergence of the series in Eq. (23) at ambient temperatures. Figure 3(a) depicts the convergence parameter *η*_{γ} as a function of the effective mass *M* for *R*_{0} = 2.8 Å and Ω = 200 cm^{−1} at *T* = 300 K. The convergence of the series in the factor of the reaction flux that originates from the quadratic term in the coupling expansion requires only two or three terms for systems with moderate to large effective force constants associated with the proton donor-acceptor mode. For systems with small effective force constants (i.e., soft hydrogen-bonded PT interfaces), full convergence may require the inclusion of more terms. For example, to achieve convergence for the 1D-HARM model with *M* = 10 amu and Ω = 200 cm^{−1} at *T* = 300 K, we had to include terms up to tenth order in *γ*. For the more realistic Morse proton potentials, the convergence parameter is sufficiently small to ensure efficient convergence even for small effective force constants.

For the thermally averaged rate constant, when an expansion of the coupling is used, the effects of the quadratic term on the rate constant can be significant for soft hydrogen-bonded PT interfaces with a small effective force constant. In this limit, the quadratic term is necessary to accurately describe the coupling at shorter proton donor-acceptor distances, which are sampled when the distribution function *P*(*R*) is broad (i.e., at high temperatures or for a very low-frequency proton donor-acceptor mode). As illustrated in Figure 3(b), however, the effects of the quadratic term become insignificant when the force constant of the proton donor-acceptor mode increases and the distribution function becomes more localized around the equilibrium value of the proton donor-acceptor distance in the reactant state. Note that the assumption of a harmonic potential along the proton donor-acceptor coordinate *R* is expected to become invalid for soft hydrogen-bonded interfaces. A more realistic anharmonic potential along the *R* coordinate would favor longer proton donor-acceptor distances, thereby avoiding the small distances at which the quadratic term in the coupling expansion is important.

For further analysis, we compared the rate constants obtained from the 2D quantum treatment of the proton donor-acceptor and proton coordinates, the thermally averaged rate constants with explicitly calculated couplings for all proton donor-acceptor distances, and the rate constants obtained from the cumulant expansion with and without the quadratic term. The results for *M* = 50 amu and *M* = 10 amu with Ω = 200 cm^{−1} are presented in Figure 4. For both models, the 2D quantum rate constants, thermally averaged rate constants with explicitly calculated couplings, and rate constants from the cumulant expansion with the quadratic term are in good agreement (dashed red, blue, and black lines, respectively). For the model with effective mass *M* = 50 amu, depicted in Figure 4(a), all of these rate constants are virtually indistinguishable for a wide temperature range. The rate constants calculated without the quadratic term in the coupling expansion are in good agreement with the others except for a slight deviation at higher temperatures (solid black line in Figure 4(a)), indicating that the effects of the quadratic term are insignificant for this model and therefore can be neglected.

The model with a lower effective mass *M* = 10 amu, corresponding to a soft hydrogen-bonded PT interface, exhibits greater differences among the various rate constant expressions, as depicted in Figure 4(b). The deviation of the rate constant calculated from the cumulant expansion with the quadratic term (dashed black line in Figure 4(b)) from the 2D quantum and thermally averaged rate constants at higher temperatures is due to the slow convergence of the cumulant expansion. These results were generated by including terms up to tenth order in *γ* in the cumulant expansion rate constant, and the inclusion of higher-order terms is expected to improve the agreement at higher temperatures. On the other hand, the slight persistent deviation of the rate constant calculated from the cumulant expansion with the quadratic term (dashed black line in Figure 4(b)) from the 2D quantum and thermally averaged rate constants at ambient and lower temperatures is due to the neglect of the dynamical interference effects in the 2D quantum and thermally averaged rate constants. In this regime, the cumulant expansion rate constant is expected to be more accurate because it includes the dynamical interference between the proton donor-acceptor mode and the energy gap.

More significant differences are observed between the rate constants calculated without the quadratic term in the coupling expansion (solid black line in Figure 4(b)) and the rate constants that include the quadratic term, implying that the quadratic term is important when shorter proton donor-acceptor distances are more accessible. Note that in realistic chemical or biological systems, the shorter proton donor-acceptor distances will not be as readily accessible even for weakly hydrogen-bonded PT interfaces due to strong repulsive nonbonding interactions between the donor and acceptor groups. Thus, the models based on a harmonic description of a low-frequency proton donor-acceptor mode appear to overemphasize the region of short distances that are not relevant for physically realistic systems.

### C. Analysis in electronically adiabatic and intermediate regimes

In the electronically adiabatic regime, typical for PT and HAT reactions, the vibronic coupling *V*(*R*) at each proton donor-acceptor distance *R* is related to the tunneling splitting for a symmetric double well proton potential on the electronically adiabatic ground state surface $U0(ad)(q,R)$. Because the models studied herein are symmetric, the proton potentials are symmetric for a fixed proton donor-acceptor distance *R*. Calculation of the semiclassical couplings and associated electron-proton adiabaticity parameters for a range of proton donor-acceptor distances for the two models LEPS-M10-Ω200 and LEPS-M50-Ω200 with different effective masses *M* = 10 amu and *M* = 50 amu, respectively, revealed that both models are in the intermediate regime. This property is illustrated in Figure 5 for the LEPS-M50-Ω200 model.

The equilibrium proton donor-acceptor distances *R*_{0} and attenuation parameters *α* and *γ* for the overlaps, tunneling splittings, and semiclassical vibronic couplings are given in Table III. The values of the attenuation parameters in these models are similar to those calculated for the overlaps in the 1D-MORSE model (Table II). Thus, the effects of the quadratic term in the coupling expansion on the rate constants in the various limits are also expected to be similar.

. | Overlap . | Tunneling splitting . | Semiclassical coupling^{a}
. | |||
---|---|---|---|---|---|---|

R_{0}, Å
. | α
. | γ
. | α
. | γ
. | α
. | γ
. |

LEPS-M10-Ω200 model | ||||||

2.72 | 17.81 | 25.43 | 16.77 | 27.67 | 16.87 | 27.23 |

LEPS-M50-Ω200 model | ||||||

2.78 | 19.75 | 25.84 | 19.33 | 26.45 | 19.39 | 26.08 |

. | Overlap . | Tunneling splitting . | Semiclassical coupling^{a}
. | |||
---|---|---|---|---|---|---|

R_{0}, Å
. | α
. | γ
. | α
. | γ
. | α
. | γ
. |

LEPS-M10-Ω200 model | ||||||

2.72 | 17.81 | 25.43 | 16.77 | 27.67 | 16.87 | 27.23 |

LEPS-M50-Ω200 model | ||||||

2.78 | 19.75 | 25.84 | 19.33 | 26.45 | 19.39 | 26.08 |

^{a}

The effective proton tunneling time is *τ*_{p} = 0.4479 (0.4267) fs, the effective electronic transition time is *τ*_{e} = 1.518 (1.518) fs, and the electron-proton adiabaticity parameter is *p* = 0.2951 (0.2811) for the LEPS-M10-Ω200 (LEPS-M50-Ω200) models.

To illustrate the effects on the thermally averaged rate constant, we examined the distance dependence of the semiclassical vibronic coupling and the approximations to this coupling based on the logarithmic expansion, as given in Eq. (4), in relation to the width of the classical distribution function for the proton donor-acceptor distances. The classical distribution function $P(R)\u221dexp\u2212\epsilon 0I(R)/kBT$ was calculated for the ground state reactant vibronic potentials $\epsilon 0I(R)$ obtained by quantization of the proton moving on the electronically diabatic reactant electronic surface *H*_{I,I}(*q*, *R*). The results are depicted in Figure 6 for two models with soft and rigid hydrogen-bonded PT interfaces. In both cases, the second-order logarithmic expansion of the coupling is sufficient to reproduce the semiclassical vibronic coupling nearly quantitatively. For the soft *R*-mode (Figure 6(b)), significant deviations of the first-order approximation to the coupling from the semiclassical vibronic coupling are observed within the sampling region covered by the distribution function *P*(*R*). In contrast, for the more rigid *R*-mode (Figure 6(a)), the first-order approximation to the vibronic coupling is satisfactory within the sampling region covered by the distribution function. Thus, the quadratic term will impact the rate constant only for the case of PT interfaces with weak hydrogen bonds and soft proton donor-acceptor vibrational modes.

As discussed above, more realistic potentials that include repulsion between the proton donor and acceptor would prevent the sampling of the unphysical short distances. In this case, the short distances sampled extensively for the soft *R*-mode system described by Figure 6(b) would be energetically inaccessible, and the probability distribution *P*(*R*) would not be a Gaussian function but rather would be steeper at short distances. As in the nonadiabatic case, the thermal averaging procedure can also be performed with the explicit calculation of the vibronic coupling for each proton donor-acceptor distance followed by numerical integration, thereby avoiding the expansion of the coupling altogether and allowing a non-Gaussian probability distribution function.

We also emphasize that the Golden Rule formalism, which is based on second-order perturbation theory, is invalid for short proton donor-acceptor distances associated with large vibronic coupling. The Golden Rule formalism is the basis for all of the rate constant expressions discussed in this paper, including those in the electronically adiabatic and nonadiabatic limits, as well as the intermediate regime. Thus, all of these nonadiabatic rate constants are expected to be invalid when the thermal averaging procedure includes significant contributions from such short proton donor-acceptor distances. Models that allow the sampling of unphysical short distances should be treated with great care within these frameworks.

## IV. CONCLUDING REMARKS

In this paper, we presented rate constant expressions for vibronically nonadiabatic PT and PCET reactions using a cumulant expansion approach based on the Golden Rule formalism. The resulting rate constant expressions are valid in both the electronically adiabatic and electronically nonadiabatic regimes, as well as the intermediate regime, as long as the appropriate form of the vibronic coupling is utilized. The general semiclassical expression for the vibronic coupling and its form in the electronically adiabatic and nonadiabatic limits are available from the literature. These rate constants differ from the previous rate constants derived with the cumulant expansion approach in that the logarithm of the vibronic coupling was expanded to second order. Thus, the logarithmic expansion of the vibronic coupling includes a quadratic as well as a linear term. The significance of this quadratic term was explored for a series of model systems.

Our analysis illustrates that inclusion of the quadratic term in the logarithmic expansion of the nonadiabatic coupling may lead to significant effects on the PT or PCET rate constant at high temperatures for PT interfaces with weak hydrogen bonds and soft proton donor-acceptor vibrational modes. The inclusion of the quadratic term in the framework of the cumulant expansion form of the nonadiabatic rate constant leads to a sign-alternating series in the expression for the time-dependent reaction flux correlation function. The convergence of this series may become slow at high temperatures for systems with a small effective force constant associated with the proton donor-acceptor mode. The effects of the quadratic term may also become significant when using the thermal averaging procedure for calculating the nonadiabatic rate constants in these regimes.

In both the cumulant expansion and the thermally averaged formulations, the effects of the quadratic term for weak hydrogen-bonding systems are less significant for more physically realistic models that prevent the sampling of unphysical short proton donor-acceptor distances. In particular, the assumption of a harmonic potential along the proton donor-acceptor coordinate is invalid for soft hydrogen-bonded interfaces, and more realistic anharmonic potentials along this coordinate inhibit sampling of the small proton donor-acceptor distances at which the quadratic term is important. Moreover, all of these rate constant expressions are based on the Golden Rule formalism, which is invalid at short proton donor-acceptor distances associated with large vibronic couplings. Thus, these rate constant expressions should be used only for systems that remain in the vibronically nonadiabatic limit (i.e., the Golden Rule regime) at all proton donor-acceptor distances sampled significantly. According to our analysis, the inclusion of the quadratic term in the coupling expansion is less likely to substantially impact the PT or PCET rate constant for physically realistic systems that remain in the Golden Rule regime.

In addition to the analysis of the quadratic term, the rigorous relation between the cumulant expansion and thermal averaging approaches was clarified. The cumulant expansion rate constant was shown to include effects arising from dynamical interference between the proton donor-acceptor and solvent motions. These dynamical effects are not included in the thermal averaging approach or in the two-dimensional full quantum approach. The general cumulant expansion and thermal averaging rate constants become mathematically equivalent for any form of the vibronic coupling when these dynamical effects are neglected. Understanding the regimes in which each rate constant expression is valid is important for future applications to PT and PCET in chemical and biological processes.

## Acknowledgments

We are grateful for support from National Institutes of Health Grant No. GM056207 (applications to enzymes) and the Center for Molecular Electrocatalysis, an Energy Frontier Research Center funded by the U.S. Department of Energy, Office of Science, Basic Energy Sciences (applications to molecular electrocatalysts).