We investigate rates of proton-coupled electron transfer (PCET) in potential sweep experiments for a generalized Anderson–Holstein model with the inclusion of a quantized proton coordinate. To model this system, we utilize a quantum classical Liouville equation embedded inside of a classical master equation, which can be solved approximately with a recently developed algorithm combining diffusional effects and surface hopping between electronic states. We find that the addition of nuclear quantum effects through the proton coordinate can yield quantitatively (but not qualitatively) different IV curves under a potential sweep compared to electron transfer (ET). Additionally, we find that kinetic isotope effects give rise to a shift in the peak potential, but not the peak current, which would allow for quantification of whether an electrochemical ET event is proton-coupled or not. These findings suggest that it will be very difficult to completely understand coupled nuclear–electronic effects in electrochemical voltammetry experiments using only IV curves, and new experimental techniques will be needed to draw inferences about the nature of electrochemical PCET.

## I. INTRODUCTION

The importance of nuclear quantum effects on electron transfer (ET) is well known in the fields of electronic gas-phase spectroscopy^{1} and electron transfer theory.^{2} For example, it is well known in gas phase spectroscopy that one must account for the displacement of nuclei when measuring electronic energies, which leads to a difference between vertical and adiabatic excitation energies.^{3} For electron transfer in solution, it is well known that a classical description of the event can be inadequate, and often, one needs to use the full quantum version of Fermi’s golden rule (FGR);^{4} there has also been an extensive discussion in the literature regarding the fact that some nuclear modes must be treated quantum mechanically, while others can be treated classically when calculating ET rates.^{5} Finally, the role of quantum nuclear effects within electrochemical ET has also long been appreciated, dating back to the original work of Butler^{6} and has been gradually expanded upon with models by Newns,^{7} Schmickler^{8} and others; here, even though the solvent may be classical, one cannot truly ignore the quantum mechanical features of fast vibrational motion on a surface. With this framework in mind, we would now like to study one more example of potential interest for strong nuclear quantum effects: the intersection of electrochemistry with proton-coupled electron transfer (PCET). There is now a wide-held belief that in many biological systems, electronic motion is quantum mechanically coupled to or gated by quantized proton motion and the same phenomena are likely to occur under electrochemical conditions.

Over the past 20 years, rate constant expressions and theories have been developed to model PCET in a wide variety of contexts.^{9–15} Yet, while much progress has been made on the theory of photoinduced PCET,^{16–20} PCET at metal oxide nanoparticle interfaces^{21} and homogeneous electrochemical PCET^{22,23} to our knowledge there have been little effort to theoretically model voltammetry experiments involving PCET. Specifically, little attention has been paid to nuclear quantum effects on electrochemical electron transfer in systems where the driving force is changed, such as in a linear sweep voltammetry experiment. In general, modeling electrochemistry (or PCET) alone is difficult enough, and so, combining classical solvent motion with ET together with quantized proton motion over long time scales with driven boundary conditions poses a substantial difficulty—both theoretical and computational.

For our part, in order to make progress, we would like to separate classical coordinates from quantum coordinates as rigorously and as generally as possible. To that end, our approach in this paper will be to invoke a partial Wigner transformation so as to map one or more quantum coordinates to classical coordinates while leaving all other coordinates unchanged. By Wigner transforming over the classical coordinates influencing the electronic motion and leaving the proton coordinate quantized, we will obtain a combined quantum classical Liouville equation-classical master equation (QCLE-CME) that can be solved with surface hopping (SH) dynamics. A similar QCLE-CME approach was used previously to study an electronic system with multiple orbitals and electronic states in the presence of a continuum of metallic electronic states.^{24} Also, recent work has shown that diffusive electrochemical ET can be modeled with a CME and simulated with SH dynamics.^{25} With this technique in mind, below, we will simulate PCET subject to a linear voltammetric sweep, which will elucidate new information about the current–potential profile of a driven PCET system.

An outline of this paper is as follows: In Sec. II, we present our model PCET Hamiltonian and show how electrochemical dynamics can be characterized by a QCLE-CME and approximately solved with SH dynamics, which we henceforth refer to as SH-CME. Because many assumptions and quite a few mathematical steps are needed to derive the relevant equations, we will be (we hope) rather pedagogical and explicit in our development of the equations. This approach will highlight how traditional ET hopping rates must be modified by Franck–Condon (i.e., overlap) factors to account for the influence of proton transfer. In Sec. III, we present results for the voltammetry simulations with different PCET parameters, illustrating how nuclear quantum effects can yield quantitatively different voltammograms. In Sec. IV, we discuss the implications of our results as far as understanding electrochemical systems that undergo PCET. In Sec. V, we conclude.

## II. THEORY AND METHODS

### A. Theoretical developments for coupled nuclear–electronic motion at a metal surface

#### 1. Anderson Holstein (AH) Hamiltonian for PCET

To model a PCET electrochemical system, imagine a redox-active molecule coupled to a bath of nuclei, whereby the molecule can exchange an electron with a metal electrode concertedly with proton transfer. We will employ a generalized Anderson–Holstein (AH) Hamiltonian of the following form:

Here, $r\u2192=[x,\zeta ,q]$ represent three relevant nuclear coordinates: the distance of the system from the electrode (*x*), a solvent reorganization coordinate (*ζ*), and a proton vibrational coordinate (*q*). The relevant momenta corresponding to those coordinates are Π_{x}, Π_{ζ}, and Π_{q}, respectively. $\psi 0el$ and $\psi 1el$ refer to electronic states of the redox active molecule, where the impurity orbital *d* is unoccupied or occupied, respectively. ${\mu}$ and ${\mu \u2032}$ refer to a complete basis of proton vibrational states for electronic states $\psi 0el$ and $\psi 1el$, respectively (for future reference). The collection of solvent nuclei is indexed by *α*, and has position and momentum *κ*_{α} and Π_{α}, respectively, with coupling to the electronic impurity given by *g*_{α}. *μ* and *μ*′ represent a given proton vibrational state when the redox molecule is uncharged and charged, respectively. *H*_{0} and *H*_{1} represent the energy of the electronic impurity orbital *d* in the unoccupied and occupied states, respectively, and the bath of metallic fermions is indexed by energy *ϵ*_{k} with coupling to the impurity orbital *d* given by *W*_{k}.

Note that the Hamiltonian in Eq. (1) is written in a slightly unorthodox way. After all, the total Hamiltonian is an operator that involves *x*, *ζ*, *q*, *k*, and *d* indices, and it acts on the tensor product vector space $Hdiff\u2297Hsolvent\u2297Hvib\u2297Hmetal\u2297Hel$. Here, $H$ signifies the Hilbert space of all functions, so, e.g., $Hel$ is spanned by the electronic impurity being occupied or not ($\psi 0$, $\psi 1$). For instance, consider the system Hamiltonian *H*_{s}, which consists of a molecule (here, just an electronic impurity level plus a proton), and a phononic bath of solvent *H*_{ph}. $U0\mu (x,\zeta )$ is a function of the diffusion and solvent coordinates so that *H*_{0} and *H*_{1} are first quantized operators that act exclusively on the vector space $Hdiff\u2297Hsolvent\u2297Hvib$. Therefore, the term $Hmol=H0\psi 0el\psi 0el+H1\psi 1el\psi 1el$ is a first quantized operator in the vector space $Hvib\u2297Hdiff\u2297Hsolvent\u2297Hel$, and so, we must assume an implicit identity operator acting on all other coordinates. At the same time, however, the term *H*_{ph} is simply written in terms of a collection of solvent coordinates and momenta. *H*_{ph} captures interactions between the nuclear motion and the solvent (as a function of distance to the electrode)—in principle, this fluctuation could also include direct interactions between the proton and solvent coordinates, although the latter are omitted below for simplicity. Finally, note that the Hamiltonian for the electronic bath ($Hbathel$) is written in second quantization in terms of some abstract electronic coupling between the system and the bath. $Hcoupe\u2212e$ is also written in second quantization, and Eq. (1h) assumes that the proton does not couple to the metallic states (for simplicity).

To reiterate, we have used a slightly unorthodox mixture of first and second quantized notation to express the Hamiltonian for a molecule that can approach a metal surface and exchange an electron (as coupled to proton motion).^{26} We have chosen this mixed notation as it will make all derivations below much simpler, as we will need to simplify the quantum Hamiltonian in Eq. (1) for any practical calculations. Nevertheless, all notation questions aside for the reader, we must now explain the key assumptions underlying the Hamiltonian in Eq. (1): *Assumption No. 1*: no direct coupling between protonic states is mediated directly by solvent [i.e., no $U0\mu \mu \u2032$ terms in Eq. (1)] and *Assumption No. 2*: no direct coupling between the proton and metallic electronic states (i.e., *ϵ*_{k} does not depend on *q*). The proton is coupled indirectly to the solvent and metal only through coupling to the molecule. Because of this choice of Hamiltonian, Eq. (1) can describe *transient* proton dynamics but not PT in its full generality for which direct solvent proton coupling must be allowed (in analogy to ET). That being said, if one sets $U0\mu =U1\mu \u2032$ (so that there is no proton reorganization), Eq. (1) clearly does include ET. More generally, however, if $U0\mu \u2260U1\mu \u2032$, Eq. (1) can describe some electronic dynamics, but all long term dynamics will necessarily include proton reorganization as well. Thus, Eq. (1) should describe concerted PCET accurately (but not sequential ET and PT).

#### 2. Redfield perturbation theory in the interaction picture

For the dynamics we are interested in studying, the quantity of interest is the reduced density matrix of the system, *ρ*_{s}, and our method for calculating *ρ*_{s} will be based on standard quantum master equation (QME) approaches.^{27,28,47,48} In coordinate form, *ρ*_{s} should depend on *q*, the impurity orbital *d*, the reaction coordinate *ζ*, and diffusion coordinate *x* so that $\rho \u2208Hmol\u2297Hmol$, where $Hmol=Hvib\u2297Hel\u2297Hdiff\u2297Hsolvent$.

Because these techniques are not standard in the *electrochemical* literature and because one must work with slightly awkward Hamiltonians, we will now present a fairly pedagogical derivation of the relevant equations. First, we begin with the exact quantum Liouville equation (QLE) for the total density matrix, *ρ*, in the interaction picture. Here, we partition the Hamiltonian such that *H*_{s} and *H*_{b} are easily diagonalized and *H*_{c} is the perturbation that is not easily diagonalizable.

Using the formal identity $\rho (t)=\u22121\u210f2\u222b0tdt\u2032\u2009Hc(t),[Hc(t\u2032),\rho (t\u2032)]$, *ρ*(*t*) can be written as

with

Second, we make the relevant QME assumptions: *Assumption No. 3a* is the assumption of weak molecule metal coupling (resulting in the Born–Markovian approximation), whereby one replaces the reduced density matrix at time *t* by $\rho (t)=\rho beq\u2297\rho s(t)$, where $\rho beq$ is the equilibrium density matrix of the bath and *ρ*_{s} is the reduced density matrix of the system. If one wants to go beyond this assumption in an *ad hoc* way, one can introduce broadening artificially, as in Ref. 29. We also assume (*Assumption No. 3b*) that the bath correlations decay fast so that $\u222b0t\u2192\u222b0\u221e$. We find, after tracing over the bath states,

At time *τ*, $Hc(\tau )=ei(Hb+Hs)\tau \u210fHc(0)e\u2212i(Hb+Hs)\tau \u210f$, and for our Hamiltonian, this expression reduces to

Note that *W*_{k} and *H*_{s} need not commute. At this point, we must unwind the double commutator in Eq. (4). [*H*_{c}(*t*), [*H*_{c}(*t* − *τ*), *ρ*(*t*)]] can be broken up into four separate commutators,

*H*_{c}(*t*−*τ*)[*H*_{c}(*t*),*ρ*(*t*)],[

*H*_{c}(*t*),*H*_{c}(*t*−*τ*)]*ρ*(*t*),−

*ρ*(*t*)[*H*_{c}(*t*),*H*_{c}(*t*−*τ*)], and−[

*H*_{c}(*t*),*ρ*(*t*)]*H*_{c}(*t*−*τ*),

yielding four distinct products

*H*_{c}(*t*)*H*_{c}(*t*−*τ*)*ρ*(*t*),−

*H*_{c}(*t*)*ρ*(*t*)*H*_{c}(*t*−*τ*),*ρ*(*t*)*H*_{c}(*t*−*τ*)*H*_{c}(*t*), and−

*H*_{c}(*t*−*τ*)*ρ*(*t*)*H*_{c}(*t*).

Noting that *H*_{c}(*τ*) can be written as the sum of two terms, this leaves eight total terms on the RHS. We start with the term No. 1a, corresponding to *H*_{c}(*t*)*H*_{c}(*t* − *τ*)*ρ*(*t*),

Here, we have used the identities $\rho s(t)=eiHst\u210f\rho s(0)e\u2212iHst\u210f$, $Trb(ck\u2020ck\rho beq)=f(\u03f5k)$ [where $f(\u03f5k)=11+e\beta \u03f5k$ is the fermi function], $\beta =1kBT$, and

Note that *ρ*_{0} + *ρ*_{1} corresponds to one particular electronic state of the molecule so that *ρ*_{0} and *ρ*_{1} are $\u2208{Hvib\u2297Hdiff\u2297Hsolvent}$. Term No. 1b corresponding to *H*_{c}(*t*)*H*_{c}(*t* − *τ*)*ρ*(*t*) becomes

Now, we will perform the same reduction on the remaining three terms. The second set of terms, corresponding to *H*_{c}(*t*)*ρ*(*t*)*H*_{c}(*t* − *τ*), are

term No. 2a:

and term No. 2b:

The third set of terms, corresponding to *ρ*(*t*)*H*_{c}(*t* − *τ*)*H*_{c}(*t*), becomes

term No. 3a:

and term No. 3b:

Finally, the fourth set of term *H*_{c}(*t* − *τ*)*ρ*(*t*)*H*_{c}(*t*) yields

term No. 4a:

term No. 4b:

This completes our derivation of the relevant equations in the interaction picture.

#### 3. Redfield perturbation theory in the Schrödinger picture

By plugging all eight terms from Eqs. (6)–(14) into Eq. (4), and converting from the interaction picture to the Schrödinger picture, we can recover two equations for the time evolution of *ρ*_{0} and *ρ*_{1},

and

These equations can be represented in the following form:

which is known as a “Redfield equation,” or equivalently

with *ρ*_{s} partitioned into *ρ*_{0} and *ρ*_{1} using Eq. (7). The superoperator $L^^$ is defined by its operation on *ρ*_{0} and *ρ*_{1},

and

where both *ρ*_{0} and *ρ*_{1} are combined electronic-vibrational density matrices.

#### 4. A partial Wigner transform

Equations (18)–(20) are still too complicated to solve numerically given that all nuclear motion is still quantum mechanical. Thus, at this point, we perform a partial Wigner transformation on Eq. (18) over the coordinates (*x*, *ζ*) while leaving the proton coordinate quantized. We will depict all quantum operators *A* with a hat, *Â*, and a superscript *vib* will signify that we have Wigner transformed over everything except the proton coordinate. Thus, we will construct $\rho 0vib(x,\zeta )$ and $\rho 1vib(x,\zeta )$. Note that $\rho 0vib$ and $\rho 1vib$ are $\u2208{Hproton}$, i.e., they are density operators in the space of the quantum proton.

When taking the Wigner transform of Eq. (18), we will need to evaluate the Wigner transform of a product of operators. Recall that $(AB)W=AWe\u210f\Lambda 2iBW$, where Λ is the Poisson bracket operator $\Lambda \u2194=\u2207\u20d6P\u22c5\u2207\u20d7R\u2212\u2207\u20d6R\u22c5\u2207\u20d7P$. Then, following the standard notation of classical mechanics, let us define ${\xd41,\xd42}=\u2211\alpha \u2202\xd41\u2202R\alpha \u2202\xd42\u2202p\alpha \u2212\u2202\xd41\u2202p\alpha \u2202\xd42\u2202R\alpha $ as the usual Poisson bracket. Note that, because *Ô*_{1} and *Ô*_{2} are operators, the order of the operators is critically important. Thus, $\u2202\xd41\u2202x\alpha \u2202\xd42\u2202p\alpha \u2260\u2202\xd41\u2202p\alpha \u2202\xd42\u2202x\alpha $. Finally, when we calculate the Wigner transformation of the commutator $\u2212i\u210f[Hs,\rho s]$ to first order in *ℏ*, which is a classical assumption about the diffusion (*x*) and solvent (*ζ*) coordinates (*Assumption No. 4*), we find the result $\u2212i\u210f[Hs,\rho s]=\u2212i\u210f[Hsvib,\rho svib]+{Hsvib,\rho svib}a$, where ${\xd41,\xd42}a=12({\xd41,\xd42}\u2212{\xd42,\xd41})$. We are left with the famous quantum classical Liouville equation (QCLE) for $\rho ^svib$,

One would like to assume that, just as for the QCLE,^{30} Eq. (21) is correct to order $O(\u210f)$. However, in Eq. (21), we have implicitly made one more (fairly standard) assumption, i.e., that $(L^^\rho s)vib=L^^vib\rho s^vib$,^{31} and so, Eq. (21) is formally less accurate than the usual QCLE. This assumption (which effectively ignores all direct nuclear quantum effects on the transient electron transfer between the metal and molecule) is largely comparable with Assumption No. 4.

Finally, at this point, we will temporarily make one final assumption, *Assumption No. 5*: we assume that *W*_{k} is independent of *q*, which will allow us to rearrange and simplify the equations for $(L^^\rho s)vib$. For a treatment of electrochemical PCET that does allow for *W*_{k}(*q*), see Ref. 32 and Appendix C. (Note that any dependence of *W*_{k} on *x*, the distance from the electrode, is different from any dependence on *q*; for a theoretical treatment of electrochemical PCET that includes the dependence of the electronic coupling on the distance *x* from the electrode, see Refs. 33 and 34.) Within the framework of Assumption No. 5, we can expand $L^^vib\rho ^0vib$ in vibrational states *μ*, *ν* [using Eqs. (19) and (20)],

where primed (unprimed) Greek indices correspond to proton vibrational states in the basis of electronic state 1 (state 0). We can similarly expand $L^^vib\rho ^1vib$ in vibrational states *μ*′, *ν*′ to obtain

At this point, if we make the substitution $\u222b0\u221ed\tau ei\u03f5k\tau \u210f=\pi \delta (\u03f5k)\u210f$, we arrive at a relatively complete expression for $\u2202\rho 0\mu \nu \u2202t$,

and similarly for $\u2202\rho 1\mu \u2032\nu \u2032\u2202t$,

Here, $\rho 0\mu \nu =\mu \rho ^0vib\nu $ and $\Gamma =2\pi \u2211kWk2\delta (\u03f5\u2212\u03f5k)$ is the metal–molecule coupling in the wide-band limit that characterizes the rate with which electrons can hop between molecule and metal.

#### 5. A secular approximation

In practice, however, it is common to make one final secular approximation, *Assumption No. 6*, whereby rapidly oscillating off-diagonal terms are neglected so that we can focus only on the diagonal (population) terms of the density matrix. *ρ*^{μμ} then becomes

with *F*_{μ→μ′} being the Franck–Condon factor (i.e., overlap integral) between wavefunctions *μ* and *μ*′, $\mu \mu \u2032$. In making the secular approximation, we ignore the quantum mechanical coherence of the proton coordinate so that we lose any capacity for protonic dynamics even at short times: we focus exclusively on coupled proton–electron transfer dynamics. Nevertheless, it will turn out that solving Eq. (26) is already a difficult enough task; to go beyond the secular approximation requires including the off-diagonal (coherence) terms, which can be difficult with surface hopping.^{35}

#### 6. Final equations

In the end, for modeling electrochemical PCET, our final equations are of the simple form

where the hopping rates are

Electrochemical PCET is clearly a generalization of electrochemical ET. If we compare Eqs. (21a) and (21b) from Ref. 25 to Eq. (28) above, the only difference is the inclusion of Franck–Condon (i.e., proton wavefunction overlap) factors in the hopping rates, and the fact that all energy differences are now between combined electron–proton vibronic energetic states, $U1\mu \u2032\u2212U0\mu $, instead of the more straightforward electronic diabatic energy difference *U*_{1} − *U*_{0}. As might be expected, hops can occur from one vibrational state (in electronic state 0 or 1) to all other possible vibrational states (in the opposite electronic state). To allow for comparison to traditional electrochemical literature definitions, we will henceforth refer to *γ*_{0→1} as *k*_{f} and *γ*_{1→0} as *k*_{b}.

Within Eqs. (27) and (28), all PCET dynamics follow directly from the potential energy surfaces $U0\mu (x,\zeta )$ and $U1\mu (x,\zeta )$, which will be recast as $UA\mu $ and $UB\mu \u2032$ (*μ* and *μ*′ index proton vibrational levels). This naming convention is to allow for more natural comparison to the chemical reaction occurring at the electrode,

where *A* represents the species oxidized (uncharged) and *B* represents the species reduced (charged). For the model calculations below, we will now make the simplifying assumption (*Assumption No. 7*) that these functions are quadratic wells over the solvent reorganization coordinate *ζ*, independent of *x*, so as to avoid a few numerical obstacles and allow for much faster calculations. In general, although, the equations can be solved with *x*-dependence. Below, we will choose our PESs to be of the form,

Here, $UA\mu $ and $UB\mu \u2032$ correspond to the vibronic contribution to the energy for electronic states A and B and proton vibrational states *μ* and *μ*′, respectively. $\u0168A\mu $ and $\u0168B\mu \u2032$ are the energy contributions for electronic states *A* and *B* (respectively) due to the quantized proton, which will depend on the nature of the proton PES and the vibrational states *μ* and *μ*′. Δ*G* is the applied driving force, and *ζ*_{A} and *ζ*_{B} are the minima positions in *ζ* for $UA\mu (\zeta )$ and $UB\mu (\zeta )$, respectively. Therefore, the solvent reorganization energy can be defined as $\lambda =12m\omega 2(\zeta B\u2212\zeta A)2$.^{36,49} In Sec. III, we will introduce different models of $\u0168$^{μ} (quadratic, Morse, and quartic double well) and demonstrate how to calculate the Franck–Condon factors for these various models. While dynamics will not be explicitly simulated over *ζ* in this work, the influence of these PESs will directly affect the hopping rates *k*_{f} and *k*_{b}, which will be discussed in Sec. II B.

### B. Simulation methods for generating IV curves on a grid

We will now describe how to use the equation developed above to model IV curves in practice.

#### 1. Classical master equation (CME) for solving the PCET AH mode in 1D

In the limit of strong friction, weak coupling Γ, and high temperature *kT* (where nuclear motions are classical), translational motion for the system is diffusive and governed by the Smoluchowski equation. Therefore, for simulating a voltammetry experiment, the dynamics from Eq. (27) can be solved by a classical master equation (CME) of the following form:

where *μ*, *μ*′ is an index for the proton vibrational state in species *A* and *B*, respectively, and $kf/b\mu \u2192\mu \u2032$ are the forward/backward hopping rates between vibrational states *μ* and *μ*′. $cA\mu $ and $cB\mu \u2032$ are the concentration densities of *A* and *B* in the *μ* and *μ*′ proton vibrational states, respectively. The *δ*(*x*) term ensures that there is hopping between electronic states only at the electrode surface at *x* = 0, and we have assumed equal diffusion coefficients for all proton vibrational states in species *A* and *B* ($DAx=DBx\u2261D$). To visualize these combined electronic-vibrational energy surfaces, see Fig. 1.

Note that Eq. (30) does not contain any *ζ* coordinate terms; we will not explicitly simulate diffusive dynamics over this coordinate, instead incorporate this dimension solely through the hopping rates $kf\mu \u2192\mu \u2032$ and $kb\mu \u2032\u2192\mu $ [see Eq. (32) below]. For all results below, we will choose spatial and initial boundary conditions as follows:

Equation (31) states that the initial concentrations at all *x* are equal to the bath concentrations of species *A* and *B*, (e) and (f) implies that these bath concentrations (at *x* = ∞) remain unchanged throughout the simulation, and (c) and (d) implies the boundary conditions at *x* = −1 are chosen to satisfy a no flux boundary condition. These boundary conditions ensure equivalence between the SH algorithm and traditional electrochemical Butler–Volmer (BV) ET boundary conditions when there are only two states involved.^{37}

Unless stated otherwise, we will employ the rate constants *k*_{f} and *k*_{b} as found by averaging the full rate constant over the solvent coordinate *ζ*, resulting in the Marcus–Hush–Chidsey (MHC) rate constant,

By reducing all dynamics to 1D and estimating solvent effects entirely through a coarse-grained hopping rate, one can achieve an enormous reduction in computational cost. For instance, we have found that simulating IV curves over both coordinates *x* and *ζ*, as was done in Ref. 25, is computationally unfeasible on our cluster with more than three proton vibrational states. For this reason, we have made *Assumption No. 7*, i.e., we assume *U*_{A} and *U*_{B} do not depend on *x*. If *U*_{A} or *U*_{B} were to depend on *x* (as well as *ζ*), a 1D reduction would not be straightforward. That being said, in Ref. 25, it was shown that simulating ET over the diffusive dimension *x* [with hopping rates given by Eq. (32)] is analytically and numerically equivalent to simulating diffusion over both *x* and *ζ* (as long as there was no explicit *x* dependence on the diabatic PESs for species *A* and *B*). Future work (on a larger computer) may well explore 2D PCET dynamics, with explicit nuclear motion in *x* and *ζ*. In Appendix A, we show that simulation in 1D with the MHC rate is equivalent to simulation in 2D for PCET in a simple case.

#### 2. Simulating voltammetry and measuring current

To simulate a voltammetry experiment, the external potential Δ*G* (which is equivalent to the difference in free energy of the product and reactant) must be changed as Eq. (30) is propagated forward in time. A common experimental approach is linear sweep voltammetry, whereby the potential is ramped linearly with scan rate *ν*,

and Δ*G*_{i} is the starting driving force. This ramp in Δ*G* is maintained until reaching the final Δ*G* value, Δ*G*_{f}. For this work, Δ*G*_{f} = −Δ*G*_{i}. To measure current, we use the difference in forward and backward rates,

where the double sum is over all proton vibrational states for both electronic states *A* and *B*. Equation (34) is equivalent to measuring the current from the flux of either species *A* or *B* at *x* = 0,

#### 3. Solving the CME: Diagonalizing the operator for diffusion and hopping

In order to propagate Eq. (30) over a long enough time to simulate a voltammetry experiment, we will diagonalize the operator for diffusion and stochastic hopping. This approach has been previously applied to study ET dynamics in a two state system with diffusion^{31} and, more recently, to study electrochemical ET under liner sweep conditions.^{25} The method below is very similar to that found in Ref. 25, with the exception being the addition of proton vibrational states through the quantized proton coordinate. In order to demonstrate the method from scratch, let us rewrite Eq. (30) in the following way:

Here, we define superoperators,

To be as specific as possible, let us now consider the case of a two vibrational state system, with energy spacing *ℏω*_{p}. Combining $LA,LB,kf$, and *k*_{b}, we can obtain the general matrix equation for *c*_{A} and *c*_{B},

The generalization to many vibrational states is straightforward. Here, in Eq. (38), the $y\u2192$ term is used to include the bath conditions from Eqs. 31(e) and 31(f), which maintain a static condition at the bath boundary, i.e., $c\u2192\u0307(x=\u221e)=0\u2192$. In other words, $y\u2192(x\u2260\u221e)=0$ and $y\u2192(x=\u221e)=\u2212(Mc\u2192)(x=\u221e)$. The size of the grid in *x* space is chosen to be large enough such that no numerical issues are introduced by the possible discontinuity of going from $y\u2192=0\u2192y\u2192=\u2212(Mc\u2192)$ at *x* = ∞.

$M$ can be non-Hermitian, since the off-diagonal elements of $M$ are not equal ($kfi\u2192i\u2032\u2260kbi\u2032\u2192i$) in general. Since diagonalization of non-Hermitian matrices is often unstable, one needs to transform $M$ into a Hermitian matrix. To do so, we define the transformation matrix *S* as

such that the new operator $M\u0303$ satisfies

Note that $M\u0303$ is Hermitian and can be easily diagonalized (also that $SLAS\u22121=LA$ because *D*_{A} is a constant, with no *x*-dependence). While the eigen*vectors* of $M\u0303$ will differ from those of $M$, $M$ and $M\u0303$ are related by a similarity transformation so that the eigen*values* are identical for both and one can analyze the dynamics of the system at any arbitrary time using either.

By diagonalizing $M\u0303$ at each time step, we can simulate a linear sweep voltammetry experiment. Once we have the eigenvectors *P* and eigenvalues Λ of $M\u0303=P\Lambda PT$, we can change variables and recast the differential equation for the time evolution of *c*_{A} and *c*_{B} into renormalized coordinates,

The final equation of motion becomes

Equation (41c) can be solved for any time step exactly given by

With Eq. (41a), we can change $b\u2192(t)$ to $c\u2192(t)=cAcB$using *P* and *S*: $c\u2192=S\u22121Pb\u2192$. Therefore, the current at each time step can be calculated with complete numerical stability.

### C. Parameters

We will work in both atomic and standard units throughout this paper and, unless stated otherwise, set *ν* = 112 V/s (10^{−16} a.u.), Γ = 0.000 27 eV (10^{−5} a.u.), *D* = 5.5 · 10^{−4} cm^{2}/s (4.75 · 10^{−4} a.u.), *ω* = 1.03 · 10^{13} s^{−1} (0.000 25 a.u.), *m* = 1.82 · 10^{−27} kg (2000 a.u.), *T* = 300 K (0.000 95 a.u.), *ζ*_{A} = −5.29 Å (−10 a.u.), *ζ*_{B} = 5.29 Å (10 a.u.), and *N*_{x} = 50. The size of the grid in *x* is given by $Lx=6DtT$, where *t*_{T} is the total time of the simulation, defined by the starting and ending Δ*G* and *ν*, $tT=|\Delta Gf\u2212\Delta Gi|\nu $, and the grid spacing is uniform such that $dx=LxNx$. In general, we will break apart a IV curve into *N*_{T} = 2000 discrete changes in driving force and use Eq. (42) to propagate the system for each discrete time step.

### D. Understanding and interpreting IV curves

Typically speaking, in electrochemistry, without any quantized proton motion, the shape of IV curves and peak parameters are functions of a unit–less parameter $\Lambda =k0D\nu FRT$.^{38,39} In the regime where Λ > 10, the kinetic behavior is reversible, i.e., the rate of ET is much faster than any other time scale in the simulation, and the peak parameters are dictated entirely by *D* and *ν*. When Λ < 0.001, the kinetic behavior is totally irreversible and the peak parameters should be dictated by the ET kinetics through the parameter *k*^{0}. Finally, when 10 ≥ Λ ≥ 0.001, the behavior is deemed “quasireversible” by Matsuda and Ayabe, where neither ET kinetics, mass transfer, or scan rate uniquely dictate the IV curve shape and peak parameters. With this framework in mind and in order to help the reader interpret our simulation data intuitively, we will include values of Λ in the caption of each figure below (for each relevant set of data). Obviously, one cannot fully understand PCET in the context of simple ET rates (see Sec. IV B), but we believe having a rough estimate of Λ should be useful as far as establishing the proper transport regime. To that end, for *k*^{0}, we use the value Γ · *FC*(0, 0), where *FC*(0, 0) is the Franck–Condon factor between ground vibrational states in electronic states 0 and 1. However, using Λ as a true estimate for the reversibility or irreversibility of a given reaction can be dangerous: after all, for the PCET framework above, for every single donor vibronic state, there are multiple acceptor vibronic states (each with a different Franck–Condon factor), and the relevant acceptor state can depend on the scan rate, diffusion constant, and/or mass (hydrogen/deuterium). Thus, one cannot ascribe to a given Hamiltonian a unique Λ value. Furthermore, we also mention that in the context of Matsuda/Yamuda Λ values, one usually presumes that electrochemical dynamics occur with Bulter–Volmer rate constants [rather than the MHC rate constants defined above in Eq. (32)]; nevertheless, we believe the Λ values reported below in each figure should be meaningful, at least qualitatively.

## III. RESULTS

All results below will depend critically on the choice of Hamiltonian for the proton, which generates the key parameters $U0\mu $ and $U1\mu \u2032$. We will work with three different functional forms of the proton potential; harmonic potentials, Morse potentials, and double well quartic potentials.

### A. Harmonic oscillator eigenstates

The simplest possible choice for proton PESs are harmonic oscillators with energy spacing *ℏω*_{p}. If *V*_{0} is the potential for the proton assuming the molecule is uncharged and *V*_{1} is the potential for the proton assuming the molecule is charged, we can write

which allows for an analytical calculation of the Franck–Condon factors. Note that $[\Pi 22mq+V0(q)]\mu =\u0168A\mu \mu $. The Franck–Condon factors are

where *ϕ*_{μ}(*q*) is the *μ*th eigenstate of the harmonic oscillator at position *q* and *g*_{p} is the normalized displacement of the two eigenstates in *q*. *F*_{μ→μ′} can be calculated by the following expression:^{40,41}

Here, p/Q are the minimum/maximum of *μ* and *μ*′, *sgn* is the sign function, and $Lmn$ is a generalized Laguerre polynomial. The reorganization energy in the proton coordinate, *λ*_{p}, is defined as $\lambda p=\u210f\omega pgp22$. In Fig. 2, we show that different IV curves do arise from inclusion of the proton coordinate. Here, (i) linear voltammograms are plotted for different combinations of *ω*_{p} and *g*_{p} along with (ii) the ratio of peak currents for PCET and ET $IpPCETIpET$ and (iii) the deviation of the peak potential between PCET and ET, $|EpPCET\u2212EpET|$. Here *E*_{p} corresponds to the overpotential at which the IV curve has its maximum value, while *I*_{p} is the maximum current attained throughout the course of the voltammetry simulation. $EpET$ and $IpET$ are the values that are obtained from a simulation with identical electronic parameters but no proton transfer. The most noticeable trends from Fig. 2 are (1) the relative independence of *I*_{p} and *E*_{p} with regard to *ω*_{p} and (2) the decrease in *I*_{p} and increase in $|EpPCET\u2212EpET|$ with regard to increase in *g*_{p}. Taken together, these two observations suggest that the main effect of the proton transfer in this model is to decrease the standard rate constant, *k*^{0}, through the incomplete orbital overlap of the two ground proton vibrational states in electronic states *A* and *B* [see Eq. (45)].

### B. Morse potential

In addition to harmonic oscillator eigenstates, the proton coordinate can be represented by a basis of Morse oscillator eigenstates. The Franck–Condon factors for these potentials can be numerically calculated using Eq. (44), with eigenstates and energy spacings from analytical solutions.^{42} In this instance, the PESs in the proton coordinate are

where $a=m\omega p22De$, and we allow *g*_{p} to vary. The energy levels are given by $En=\u210f\omega p2De+4nDe\u2212\u210f\omega p(n2+n+14)4De$, and we set *D*_{e} = 10*ω*_{p}. Therefore, *ω*_{p} affects the energy level spacing, and *g*_{p} affects the proton wavefunction overlap between states *A* and *B*. The voltammograms comparing ET to PCET with Morse oscillators for the proton states are shown in Fig. 4. We see some of the same qualitative behavior as compared with the harmonic oscillator case, whereby *I*_{p} decreases and *E*_{p} increases as a function of increased spacing in the proton coordinate. However, Fig. 4(d) shows that at large spacings in *q*, the peak potential is dependent on *ω*_{p}, in contrast to the harmonic oscillator case. This suggests that the ET mechanism can be different in the Morse oscillator case vs the harmonic oscillator case; apparently, hops can occur (in practice) between the ground vibrational state in *A* and excited vibrational states in *B*.

### C. Quartic potential

Our last model potential will be the quartic double well proton PES.^{43–45} The PESs are given by

where *q*_{c} is the crossing of the two PESs for *A* and *B*. For our simulations, we restrict ourselves to the case where the *a* coefficients are the same for species *A* and *B* (i.e., $aiA=aiB$ for *i* = 1, 2, 4). *L* is a reference length scale, $L=\u210fm\omega p$, and *E*_{ref} is a reference energy scale. *q*_{c} is chosen as $a2\u210f4a4m\omega p$, and shown in Fig. 5 are example quartic PESs for four different values of *a*_{1}. When *a*_{4} = *a*_{1} = 0, *a*_{2} = −1, this model is identical to the quadratic harmonic oscillator model. The parameter *a*_{1} dictates the asymmetry of the two wells in the quartic model and yields symmetric double wells when *a*_{1} = 0. We compute the energy levels/spacings and wavefunctions/Franck–Condon factors numerically, and we use the lowest ten levels in our simulations. In Fig. 6, we show voltammograms for different values of *ω*_{p} and *a*_{1}. When *a*_{1} = 0, the proton coordinate PESs are identical for both electronic states *A* and *B*, and as expected, we see no deviation from the ET case for this parameter choice. However, as *a*_{1} is increased, the asymmetry between proton PESs increases, introducing quantitative differences in *I*_{p} and *E*_{p} when compared to ET. Similar to the other two models examined, *I*_{p} decreases and *E*_{p} increases as a function of increasing *g*_{p} (i.e., decreasing the orbital overlap between ground proton vibrational states). However, much like the Morse potential, *I*_{p} and *E*_{p} both depend on *ω*_{p}, which is also a direct probe of the energy spacings in this model.^{46} Interestingly, these results, along with those from the Morse potential, suggest that including anharmonicity in the proton PES can result in different mechanisms for ET when including a quantized proton.

### D. Proton coordinate spacing of *A* and *B* as an additional degree of freedom

Up to this point, we have restricted our attention to cases where the displacement between the proton PESs $\u0168A\mu $ and $\u0168B\mu \u2032$ is constant (see Fig. 2). However, it is well known from the PCET literature that this spacing in *q* between energetic minima can change on a very fast timescale due to the proton donor–acceptor vibrational motion, which can be influenced by other coordinates within the system as well. One way to model a non-constant displacement in the proton PESs is by sampling from a distribution of displacement values at each time step. To this end, we can perform simulations in the harmonic oscillator basis where *g*_{p} [see Eq. (44)] at every time step of the simulation is drawn from a lognormal distribution of *g*_{p} values, with mean $g\xaf$_{p} and variance *σ*. Plotted in Fig. 7 are voltammograms for various *g*_{p} distribution parameters, along with the same simulations for a single *g*_{p} value, equal to $g\xaf$_{p}. Overall, the results have a simple interpretation. At small *σ* values, the distribution of *g*_{p} values simulations gives the same IV curve as the single *g*_{p} value, as anticipated. However, at larger *σ* values, we see an earlier onset of current with respect to overpotential. This early onset effect is due to the presence of smaller *g*_{p} values in our electrochemical simulations, allowing for the possibility of much larger Franck–Condon factors between the ground vibrational states in *A* and *B* as compared to the case when *g*_{p} is fixed. This finding is particularly noteworthy, as experimental IV curves are more likely to resemble the noisy curves seen in Fig. 7 than those of the static *g*_{p} case in Fig. 2.

## IV. DISCUSSION

The data above have demonstrated that the presence of proton motion can influence the nature of the IV curve, leading to final results that can be quantitatively different from the bare ET case. For instance, we observe that *E*_{p} can often occur at large overpotentials that depend on the nature of the proton PES and *I*_{p} can be reduced by as much as 40%. These differences between ET and PCET behavior represent a significant finding. As such, one would very much like a diagnostic tool for assessing whether or not proton motion is inherent in a given electrochemical experiment. With this intention in mind, one can ask whether one can find single quantitative diagnostics of PCET. To that end, we have investigated the usual electrochemical tools from the literature: first, the dependence of peak current and peak potential on diffusion strength (as measured by the diffusion constant) and scan rate and second, the effect on the IV curve of changing the hydrogen isotope to deuterium [related to the kinetic isotope effect (KIE)]. While we find that scan rate and diffusion constant analysis do not yield conclusive evidence of PCET, it appears that isotope effects have a specific signature, whereby the peak potential *E*_{p} can shift as a function of isotope but *I*_{p} remains relatively unchanged.

### A. Dependence on diffusion constant and scan rate

In Fig. 8, we show the dependence of peak current and peak potential on diffusion constant and scan rate for multiple *g*_{p} values (in the harmonic oscillator eigenstate PCET model). The *g*_{p} = 0 trace is the same as the ET case and is included for reference. We see that the peak current closely follows the expected square root dependence on the diffusion constant and scan rate for all values of *g*_{p} studied.^{39} For *E*_{p}, there appears to be a slight change in the behavior with respect to *g*_{p}; however, because there is no simple analytical form for *E*_{p} that we can compare against, it is not clear whether plotting *E*_{p} vs scan rate or diffusion constant can serve as a useful means of identifying whether a quantized proton is important for the mechanism being investigated.

### B. Hydrogen/deuterium isotope effects in a quartic potential

Let us now discuss another important effect in PCET, the kinetic isotope effect (KIE), whereby different isotopes of the hydrogen being transferred can greatly influence the rate constant of the overall reaction. The effect of changing the isotope from hydrogen to deuterium leads to different Franck–Condon factors and effective coupling elements and can lead to large changes in rate constants in other contexts. As far as electrochemical PCET simulations are concerned, within our quartic model, changing the isotope of hydrogen amounts to simply changing the mass and thereby the energy levels and the protonic wavefunctions.

The first thing we note is that, apparently, the peak current *I*_{p} does not seem to differentiate much between hydrogen and deuterium. This creates a very interesting state of affairs. After all, in Fig. 2 (for the case of harmonic oscillators), we showed that *I*_{p} can depend on *g*_{p} (and the electron transfer rate) in the reversible or quasi-reversible regimes. In Fig. 9, however, we show that *I*_{p} can sometime be roughly independent of *g*_{p} in this same regime—the relative contribution of donor and acceptor states may change as a function of mass (just as they change as a function of *E*_{ref} in Fig. 6). Obviously, one cannot perfectly match up intuition from the standard electrochemical theory (with only two quantized, electronic states) with the current problems in PCET (with more than two quantized, vibronic states): the notion of reversibility and irreversibility becomes more complicated (also interesting). To complicate the picture even more, note that we are using Marcus–Hush rate constants, rather than the more standard Butler–Volmer rate constants (see Sec. II B 1).

Finally, let us turn to the peak potential. In Fig. 9, we investigate IV curves that arise when we change the proton to deuterium, and we plot voltammograms for multiple *a*_{1} values from the quartic potential model, as well as a plot of $EpPCET(D)\u2212EpPCET(H)$. For the peak potential, hydrogen exhibits a much earlier peak current than deuterium on the order of up to 150 mV. Thus, one possible signature of PCET would be an IV curve whose peak position is shifted but the peak current is not reduced significantly. While these results on their own do not yet offer any rigorous means of quantifying the absolute effect of isotope on the IV curve for PCET vs ET, they do suggest that the identity of the isotope being transferred concertedly with ET can have a measurable effect on the resultant voltammogram.

## V. CONCLUSION

In this paper, we have presented a means of simulating PCET dynamics under a time dependent driving force by propagating a SH-CME, which can be solved with SH dynamics. Our approach is based on the following seven assumptions:

Assumption No. 1: no direct coupling between the proton and the solvent [no $U0\mu \mu \u2032$ terms in Eq. (1)].

Assumption No. 2: no direct coupling between the proton and metallic electronic states (

*ϵ*_{k}does not depend on*q*).Assumption No. 3: Born–Markovian approximation (weak metal–molecule coupling and fast metallic electrons); $\rho (t)=\rho beq\u2297\rho s(t)$.

Assumption No. 4: Classical

*x*and*ζ*.Assumption No. 5:

*W*_{k}(electronic coupling) is independent of*q*(proton coordinate).Assumption No. 6: Secular approximation (off-diagonal terms of density matrix ignored).

Assumption No. 7:

*U*_{A}and*U*_{B}(electronic PESs) are independent of*x*(distance from the electrode).

We find that merging PT and ET in this system yields quantitatively different voltammograms, although the IV curves are often similar qualitatively (just shifted); the exact resultant effect depends on the specific nature of the proton PES or excited vibrational states, and from the data above, we find it can be very difficult to know ahead of time whether one is operating in the reversible or irreversible regimes. Clearly, with more than one acceptor and more than one donor vibronic state, the entire concept of reversible and irreversible regimes becomes more complicated. New analytic theory will be needed to explore the nature of *I*_{p}/*E*_{p} curves as a function of a collection of reduced parameters along the lines of what Matsuda did for standard electrochemical dynamics. Finally, turning to the specific question of isotopes, we find that isotope effects lead to deviations in *E*_{p} (although not *I*_{p}). Ultimately, however, we find that an IV curve is most likely insufficient on its own to establish the presence of PCET and will need to be supplemented by some other type of experiment for characterizing the nature of electrochemical dynamics.

Now, some of the observations above are dependent on our choice of model Hamiltonian; namely, the lack of direct interaction between the proton and the solvent and an absence of any relaxation pathways between excited vibrational states on a given electronic state. This choice of Hamiltonian thus prevents us from seeing alternative dynamics, such as sequential proton and electron transfer. Looking forward, the next step would be to perform these simulations with actual PCET PESs inferred from theory and experiment; once this is done, parameters could be obtained from *ab initio* calculations and this framework could be applied to known experimental PCET systems to see if this approach can recover experimentally observed IV curves that deviate from current electrochemical simulation models. Also, pH effects can have a large influence on the dynamics of electrochemical systems undergoing PCET events, which were neglected for this work. Looking forward, several arenas still remain unexplored in our quest to understand PCET in electrochemical dynamics. For the interested reader, the code used in this paper will soon be available on GitHub at the address https://github.com/alecja/electrochemistry_sim.

## ACKNOWLEDGMENTS

We thank the DoD High Performance Computing Modernization Program for computer time. This work was supported by the U.S. Air Force Office of Scientific Research (USAFOSR) under Grant Nos. FA9550-18-1-0497 and FA9550-18-1-0420.

### APPENDIX A: EQUIVALENCE OF 2D AND 1D SIMULATIONS

In Ref. 25, it was shown that, in the absence of PES entanglement between coordinates, a 1D model with Marcus–Hush–Chidsey rates gave quantitatively identical results to simulations over 2D (with a solvent coordinate explicitly modeled in addition to a diffusive coordinate). We would like to confirm that this result still stands true for PCET and that inclusion of the proton coordinate does not break this equivalency between 1D and 2D models. In Fig. 10, we show the equivalence of full 2D simulations over *x* and *ζ* with 1D simulations over only *x* for the harmonic oscillator model, where the contribution of *ζ* is incorporated through the Marcus–Hush–Chidsey rate. We see quantitative agreement for all values of *ω*_{p}, suggesting that the 1D model [in Eq. (30)] is sufficient for capturing the IV behavior in these simulations.

### APPENDIX B: CONVERGENCE OF SIMULATIONS WITH RESPECT TO NUMBER OF VIBRATIONAL STATES

Throughout this paper, we have limited our attention to *N*_{μ} = 2 for the harmonic oscillator system, but one may wonder whether different effects would emerge if we were to include even more vibrational states. Plotted in Fig. 11 are simulations for *N*_{μ} = 2 and *N*_{μ} = 5 for multiple combinations of *ω*_{p} and *g*_{p}. Since the two simulations give identical results, we can surmise that there is no additional information lost or gained by choosing *N*_{μ} anywhere in the range from 2 to 5, and all simulations seem converged by *N*_{μ} = 5. This conclusion makes intuitive sense, after all, as for even the smallest vibrational state spacing examined (*ω*_{p} = 8.42*kT*), there is a spacing of over 42*kT* between the ground and highest excited vibrational state.

### APPENDIX C: *W*_{k} DEPENDENCE ON PROTON COORDINATE

The QCLE-CME derivation of Sec. II can be expanded to allow for dependence of the molecule–metal coupling, *W*_{k}, on the proton coordinate *q*. The following derivation will focus on *ρ*_{0} only, as the final result for *ρ*_{1} can easily be solved by analogy. We first start with the analog Eq. (22),

Next, we arrive at the analog of Eq. (24),

Now, we define a wide-band limit molecule–metal coupling strength that is dependent on the proton vibrational state,

which when is substituted into Eq. (C2) yields

Finally, taking the secular approximation, we arrive at

which is similar to Eq. (26), except that $\Gamma F\mu \u2192\mu \u20322$ has been replaced with the more general $\Gamma \mu \u2032\mu \mu \mu \u2032$.

## REFERENCES

*H*_{sys} can also be expressed in second quantized notation as $Hsys=\u2211\mu U0\mu (1\u2212d\u2020d)a\mu \u2020a\mu +\u2211\mu \u2032U1\mu \u2032d\u2020d\u2009a\mu \u2032\u2020a\mu \u2032+Hph$, where $a\mu (a\mu \u2020)$ are annihilation (creation) operators for the proton vibrational state *μ* and $1$ is the identity operator.

Of course, *W*_{k} must depend strongly on *x*, and it is for this reason that we allow hopping only at the *x* = 0 grid point. Moreover, in Ref. 25, it was shown that, as long as the coupling decays quickly with distance from the electrode (e.g. exponentially), we can incorporate an x-dependence in the coupling by integrating over the relevant hopping regime and redefining the hopping rate at the first grid point in *x*. In practice, for any reasonably spaced uniform grid (which is usually on the order of microns), all relevant hopping should be only at the first grid point.

We note that recent work by Matyushov and Newton^{49} has shown that nonergodic effects due to a slowly relaxing medium (of comparable order to the reaction time) can play a significant role in electrochemical ET by reducing the solvent reorganization energy below its thermodynamic limit from Marcus theory. This effect is mainly noticeable in electrochemical ET due to the multiple orders of magnitude in rates covered over the course of the entire potential sweep. While we do not account for this effect in our framework, this observation provides an additional way that one could fine tune the above model by incorporating an overpotential dependence into λ.

For a more rigorous justification for this choice of electrode boundary condition, see Ref. 25.

*a*

_{4}and

*a*

_{2}can also change the observed IV behavior but do not have as simple of physical intuition as the

*a*

_{1}parameters with regards to the shape of the proton PESs for electronic states A and B. Thus, we have chosen to show only the

*a*

_{1}dependence, as this illustrates qualitatively how the quartic oscillator shows similar current vs potential trends as the other oscillators examined.