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.
The importance of nuclear quantum effects on electron transfer (ET) is well known in the fields of electronic gas-phase spectroscopy1 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 Butler6 and has been gradually expanded upon with models by Newns,7 Schmickler8 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 interfaces21 and homogeneous electrochemical PCET22,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, 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. and refer to electronic states of the redox active molecule, where the impurity orbital d is unoccupied or occupied, respectively. and refer to a complete basis of proton vibrational states for electronic states and , 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. H0 and H1 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 Wk.
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 . Here, signifies the Hilbert space of all functions, so, e.g., is spanned by the electronic impurity being occupied or not (, ). For instance, consider the system Hamiltonian Hs, which consists of a molecule (here, just an electronic impurity level plus a proton), and a phononic bath of solvent Hph. is a function of the diffusion and solvent coordinates so that H0 and H1 are first quantized operators that act exclusively on the vector space . Therefore, the term is a first quantized operator in the vector space , and so, we must assume an implicit identity operator acting on all other coordinates. At the same time, however, the term Hph is simply written in terms of a collection of solvent coordinates and momenta. Hph 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 () is written in second quantization in terms of some abstract electronic coupling between the system and the bath. 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 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 (so that there is no proton reorganization), Eq. (1) clearly does include ET. More generally, however, if , 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 , where .
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 Hs and Hb are easily diagonalized and Hc is the perturbation that is not easily diagonalizable.
Using the formal identity , ρ(t) can be written as
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 , where 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 . We find, after tracing over the bath states,
At time τ, , and for our Hamiltonian, this expression reduces to
Note that Wk and Hs need not commute. At this point, we must unwind the double commutator in Eq. (4). [Hc(t), [Hc(t − τ), ρ(t)]] can be broken up into four separate commutators,
Hc(t − τ)[Hc(t), ρ(t)],
[Hc(t), Hc(t − τ)]ρ(t),
−ρ(t)[Hc(t), Hc(t − τ)], and
−[Hc(t), ρ(t)]Hc(t − τ),
yielding four distinct products
Hc(t)Hc(t − τ)ρ(t),
−Hc(t)ρ(t)Hc(t − τ),
ρ(t)Hc(t − τ)Hc(t), and
−Hc(t − τ)ρ(t)Hc(t).
Noting that Hc(τ) 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 Hc(t)Hc(t − τ)ρ(t),
Here, we have used the identities , [where is the fermi function], , and
Note that ρ0 + ρ1 corresponds to one particular electronic state of the molecule so that ρ0 and ρ1 are . Term No. 1b corresponding to Hc(t)Hc(t − τ)ρ(t) becomes
Now, we will perform the same reduction on the remaining three terms. The second set of terms, corresponding to Hc(t)ρ(t)Hc(t − τ), are
term No. 2a:
and term No. 2b:
The third set of terms, corresponding to ρ(t)Hc(t − τ)Hc(t), becomes
term No. 3a:
and term No. 3b:
Finally, the fourth set of term Hc(t − τ)ρ(t)Hc(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
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 is defined by its operation on ρ0 and ρ1,
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 and . Note that and are , 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 , where Λ is the Poisson bracket operator . Then, following the standard notation of classical mechanics, let us define as the usual Poisson bracket. Note that, because Ô1 and Ô2 are operators, the order of the operators is critically important. Thus, . Finally, when we calculate the Wigner transformation of the commutator to first order in ℏ, which is a classical assumption about the diffusion (x) and solvent (ζ) coordinates (Assumption No. 4), we find the result , where . We are left with the famous quantum classical Liouville equation (QCLE) for ,
One would like to assume that, just as for the QCLE,30 Eq. (21) is correct to order . However, in Eq. (21), we have implicitly made one more (fairly standard) assumption, i.e., that ,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 Wk is independent of q, which will allow us to rearrange and simplify the equations for . For a treatment of electrochemical PCET that does allow for Wk(q), see Ref. 32 and Appendix C. (Note that any dependence of Wk 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 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 in vibrational states μ′, ν′ to obtain
At this point, if we make the substitution , we arrive at a relatively complete expression for ,
and similarly for ,
Here, and 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 μ′, . 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, , instead of the more straightforward electronic diabatic energy difference U1 − U0. 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 kf and γ1→0 as kb.
Within Eqs. (27) and (28), all PCET dynamics follow directly from the potential energy surfaces and , which will be recast as and (μ 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, and correspond to the vibronic contribution to the energy for electronic states A and B and proton vibrational states μ and μ′, respectively. and 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 and , respectively. Therefore, the solvent reorganization energy can be defined as .36,49 In Sec. III, we will introduce different models of μ (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 kf and kb, 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 are the forward/backward hopping rates between vibrational states μ and μ′. and 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 (). 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 and [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 kf and kb 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 UA and UB do not depend on x. If UA or UB 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 ΔGi is the starting driving force. This ramp in ΔG is maintained until reaching the final ΔG value, ΔGf. For this work, ΔGf = −ΔGi. 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 diffusion31 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 , and kb, we can obtain the general matrix equation for cA and cB,
The generalization to many vibrational states is straightforward. Here, in Eq. (38), the 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., . In other words, and . 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 at x = ∞.
can be non-Hermitian, since the off-diagonal elements of are not equal () in general. Since diagonalization of non-Hermitian matrices is often unstable, one needs to transform into a Hermitian matrix. To do so, we define the transformation matrix S as
such that the new operator satisfies
Note that is Hermitian and can be easily diagonalized (also that because DA is a constant, with no x-dependence). While the eigenvectors of will differ from those of , and are related by a similarity transformation so that the eigenvalues are identical for both and one can analyze the dynamics of the system at any arbitrary time using either.
By diagonalizing at each time step, we can simulate a linear sweep voltammetry experiment. Once we have the eigenvectors P and eigenvalues Λ of , we can change variables and recast the differential equation for the time evolution of cA and cB 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 to using P and S: . Therefore, the current at each time step can be calculated with complete numerical stability.
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 cm2/s (4.75 · 10−4 a.u.), ω = 1.03 · 1013 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 Nx = 50. The size of the grid in x is given by , where tT is the total time of the simulation, defined by the starting and ending ΔG and ν, , and the grid spacing is uniform such that . In general, we will break apart a IV curve into NT = 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 .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 k0. 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 k0, 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.
All results below will depend critically on the choice of Hamiltonian for the proton, which generates the key parameters and . 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 V0 is the potential for the proton assuming the molecule is uncharged and V1 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 . The Franck–Condon factors are
where ϕμ(q) is the μth eigenstate of the harmonic oscillator at position q and gp 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 is a generalized Laguerre polynomial. The reorganization energy in the proton coordinate, λp, is defined as . 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 gp along with (ii) the ratio of peak currents for PCET and ET and (iii) the deviation of the peak potential between PCET and ET, . Here Ep corresponds to the overpotential at which the IV curve has its maximum value, while Ip is the maximum current attained throughout the course of the voltammetry simulation. and 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 Ip and Ep with regard to ωp and (2) the decrease in Ip and increase in with regard to increase in gp. Taken together, these two observations suggest that the main effect of the proton transfer in this model is to decrease the standard rate constant, k0, 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 , and we allow gp to vary. The energy levels are given by , and we set De = 10ωp. Therefore, ωp affects the energy level spacing, and gp 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 Ip decreases and Ep 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 qc 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., for i = 1, 2, 4). L is a reference length scale, , and Eref is a reference energy scale. qc is chosen as , and shown in Fig. 5 are example quartic PESs for four different values of a1. When a4 = a1 = 0, a2 = −1, this model is identical to the quadratic harmonic oscillator model. The parameter a1 dictates the asymmetry of the two wells in the quartic model and yields symmetric double wells when a1 = 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 a1. When a1 = 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 a1 is increased, the asymmetry between proton PESs increases, introducing quantitative differences in Ip and Ep when compared to ET. Similar to the other two models examined, Ip decreases and Ep increases as a function of increasing gp (i.e., decreasing the orbital overlap between ground proton vibrational states). However, much like the Morse potential, Ip and Ep 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 and 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 gp [see Eq. (44)] at every time step of the simulation is drawn from a lognormal distribution of gp values, with mean p and variance σ. Plotted in Fig. 7 are voltammograms for various gp distribution parameters, along with the same simulations for a single gp value, equal to p. Overall, the results have a simple interpretation. At small σ values, the distribution of gp values simulations gives the same IV curve as the single gp 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 gp 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 gp 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 gp case in Fig. 2.
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 Ep can often occur at large overpotentials that depend on the nature of the proton PES and Ip 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 Ep can shift as a function of isotope but Ip 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 gp values (in the harmonic oscillator eigenstate PCET model). The gp = 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 gp studied.39 For Ep, there appears to be a slight change in the behavior with respect to gp; however, because there is no simple analytical form for Ep that we can compare against, it is not clear whether plotting Ep 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 Ip 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 Ip can depend on gp (and the electron transfer rate) in the reversible or quasi-reversible regimes. In Fig. 9, however, we show that Ip can sometime be roughly independent of gp 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 Eref 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 a1 values from the quartic potential model, as well as a plot of . 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.
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 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); .
Assumption No. 4: Classical x and ζ.
Assumption No. 5: Wk (electronic coupling) is independent of q (proton coordinate).
Assumption No. 6: Secular approximation (off-diagonal terms of density matrix ignored).
Assumption No. 7: UA and UB (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 Ip/Ep 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 Ep (although not Ip). 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.
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 gp. 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.42kT), there is a spacing of over 42kT between the ground and highest excited vibrational state.
APPENDIX C: Wk DEPENDENCE ON PROTON COORDINATE
The QCLE-CME derivation of Sec. II can be expanded to allow for dependence of the molecule–metal coupling, Wk, 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 has been replaced with the more general .
Hsys can also be expressed in second quantized notation as , where are annihilation (creation) operators for the proton vibrational state μ and is the identity operator.
Of course, Wk 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 Newton49 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 λ.