Modeling voltammetry curves for proton coupled electron transfer: The importance of nuclear quantum effects

Coffman, Alec J.; Dou, Wenjie; Hammes-Schiffer, Sharon; Subotnik, Joseph E.

doi:10.1063/5.0010412

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¹ and electron transfer theory.² 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.³ 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);⁴ 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.⁵ Finally, the role of quantum nuclear effects within electrochemical ET has also long been appreciated, dating back to the original work of Butler⁶ and has been gradually expanded upon with models by Newns,⁷ Schmickler⁸ 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²¹ 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.²⁴ Also, recent work has shown that diffusive electrochemical ET can be modeled with a CME and simulated with SH dynamics.²⁵ 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:

H = H_{s} + H_{b} + H_{c},

(1a)

H_{s} = H_{s y s} \equiv H_{m o l} + H_{p h},

(1b)

H_{s y s} = \underset{H_{m o l}}{\underset{︸}{H_{0} |ψ_{0}^{e l}⟩ ⟨ψ_{0}^{e l}| + H_{1} |ψ_{1}^{e l}⟩ ⟨ψ_{1}^{e l}|}} + H_{p h},

(1c)

H_{0} = \sum_{μ} U_{0}^{μ} (x, ζ) |μ⟩ ⟨μ|,

(1d)

H_{1} = \sum_{μ^{'}} U_{1}^{μ^{'}} (x, ζ) |μ^{'}⟩ ⟨μ^{'}|,

(1e)

\begin{align} H_{p h} = \sum_{α} [\frac{1}{2} (\frac{Π_{α}^{2}}{m_{α}} + m_{α} ω_{α}^{2} (κ_{α} - \frac{g_{α}^{x} x}{m_{α} ω_{α}^{2}} - \frac{g_{α}^{ζ} ζ}{m_{α} ω_{α}^{2}}))] + \frac{Π_{ζ}^{2}}{2 m_{ζ}} + \frac{Π_{x}^{2}}{2 m_{x}}, \end{align}

(1f)

H_{b} = H_{b a t h}^{e l} = \sum_{k} ϵ_{k} c_{k}^{†} c_{k},

(1g)

H_{c} = H_{c o u p}^{e - e} = \sum_{k} W_{k} (\vec{r}) (c_{k}^{†} d + d^{†} c_{k}) .

(1h)

Here, $\vec{r} = [x, ζ, 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. $|ψ_{0}^{e l}⟩$ and $|ψ_{1}^{e l}⟩$ 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 $ψ_{0}^{e l}$ and $ψ_{1}^{e l}$ ⁠, 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₀ and H₁ 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 $H_{d i f f} \otimes H_{s o l v e n t} \otimes H_{v i b} \otimes H_{m e t a l} \otimes H_{e l}$ ⁠. Here, $H$ signifies the Hilbert space of all functions, so, e.g., $H_{e l}$ is spanned by the electronic impurity being occupied or not (⁠ $|ψ_{0}⟩$ ⁠, $|ψ_{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. $U_{0}^{μ} (x, ζ)$ is a function of the diffusion and solvent coordinates so that H₀ and H₁ are first quantized operators that act exclusively on the vector space $H_{d i f f} \otimes H_{s o l v e n t} \otimes H_{v i b}$ ⁠. Therefore, the term $H_{m o l} = H_{0} |ψ_{0}^{e l}⟩ ⟨ψ_{0}^{e l}| + H_{1} |ψ_{1}^{e l}⟩ ⟨ψ_{1}^{e l}|$ is a first quantized operator in the vector space $H_{v i b} \otimes H_{d i f f} \otimes H_{s o l v e n t} \otimes H_{e l}$ ⁠, 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 (⁠ $H_{b a t h}^{e l}$ ⁠) is written in second quantization in terms of some abstract electronic coupling between the system and the bath. $H_{c o u p}^{e - e}$ 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).²⁶ 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 $U_{0}^{μ μ^{'}}$ 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 $U_{0}^{μ} = U_{1}^{μ^{'}}$ (so that there is no proton reorganization), Eq. (1) clearly does include ET. More generally, however, if $U_{0}^{μ} \neq U_{1}^{μ^{'}}$ ⁠, 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 $ρ \in H_{m o l} \otimes H_{m o l}$ ⁠, where $H_{m o l} = H_{v i b} \otimes H_{e l} \otimes H_{d i f f} \otimes H_{s o l v e n t}$ ⁠.

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 $ρ (t) = - \frac{1}{ℏ^{2}} \int_{0}^{t} d t^{'} [) H_{c} (t), [H_{c} (t^{'}), ρ (t^{'})]$ ⁠, ρ(t) can be written as

\frac{\partial ρ (t)}{\partial t} = - \frac{i}{ℏ} [H_{c}, ρ (0)] - \frac{1}{ℏ^{2}} \int_{0}^{t} d t^{'} [H_{c} (t), [H_{c} (t^{'}), ρ (t^{'})]],

(2)

with

H_{c} (t) \equiv e^{\frac{i (H_{s} + H_{b}) t}{ℏ}} H_{c} e^{\frac{- i (H_{s} + H_{b}) t}{ℏ}} .

(3)

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 $ρ (t) = ρ_{b}^{e q} \otimes ρ_{s} (t)$ ⁠, where $ρ_{b}^{e q}$ 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 $\int_{0}^{t} \to \int_{0}^{\infty}$ ⁠. We find, after tracing over the bath states,

\frac{\partial ρ_{s} (t)}{\partial t} = - \frac{1}{ℏ^{2}} \int_{0}^{\infty} d τ T r_{b} [H_{c} (t), [H_{c} (t - τ), ρ_{b}^{e q} \otimes ρ_{s} (t)]] .

(4)

At time τ, $H_{c} (τ) = e^{\frac{i (H_{b} + H_{s}) τ}{ℏ}} H_{c} (0) e^{\frac{- i (H_{b} + H_{s}) τ}{ℏ}}$ ⁠, and for our Hamiltonian, this expression reduces to

\begin{align} H_{c} (τ) = & \sum_{k} e^{\frac{i (H_{b} + H_{s}) τ}{ℏ}} W_{k} (\vec{r}) (c_{k}^{†} d + d^{†} c_{k}) e^{- \frac{i (H_{b} + H_{s}) τ}{ℏ}} \\ = & \sum_{k} e^{\frac{i ϵ_{k} τ}{ℏ}} c_{k}^{†} e^{\frac{i H_{s} τ}{ℏ}} W_{k} (\vec{r}) d e^{- \frac{i H_{s} τ}{ℏ}} + e^{\frac{i H_{s} τ}{ℏ}} W_{k} (\vec{r}) d^{†} e^{- \frac{i H_{s} τ}{ℏ}} c_{k} e^{\frac{- i ϵ_{k} τ}{ℏ}} . \end{align}

(5)

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),

\begin{align} T r_{b} \sum_{k, k^{'}} (e^{\frac{i ϵ_{k} t}{ℏ}} e^{- \frac{i ϵ_{k^{'}} (t - τ)}{ℏ}} c_{k}^{†} c_{k}^{'} e^{\frac{i H_{s} t}{ℏ}} W_{k} (\vec{r}) d e^{- \frac{i H_{s} t}{ℏ}}) \\ (\times e^{\frac{i H_{s} (t - τ)}{ℏ}} W_{k^{'}} (\vec{r}) d^{†} e^{- \frac{i H_{s} (t - τ)}{ℏ}}) ρ_{b}^{e q} \otimes ρ_{s} (t) \\ = \sum_{k} e^{\frac{i ϵ_{k} τ}{ℏ}} f (ϵ_{k}) e^{\frac{i H_{s} t}{ℏ}} W_{k} (\vec{r}) d e^{- \frac{i H_{s} τ}{ℏ}} W_{k} (\vec{r}) d^{†} e^{- \frac{i H_{s} (t - τ)}{ℏ}} \\ \times e^{\frac{i H_{s} t}{ℏ}} ρ_{s} (0) e^{- \frac{i H_{s} t}{ℏ}} \\ = \sum_{k} e^{\frac{i ϵ_{k} τ}{ℏ}} f (ϵ_{k}) e^{\frac{i H_{s} t}{ℏ}} W_{k} (\vec{r}) e^{- \frac{i H_{1} τ}{ℏ}} W_{k} (\vec{r}) e^{\frac{i H_{0} τ}{ℏ}} ρ_{0} |ψ_{0}^{e l}⟩ ⟨ψ_{0}^{e l}| e^{- \frac{i H_{s} t}{ℏ}} . \end{align}

(6)

Here, we have used the identities $ρ_{s} (t) = e^{\frac{i H_{s} t}{ℏ}} ρ_{s} (0) e^{- \frac{i H_{s} t}{ℏ}}$ ⁠, $T r_{b} (c_{k}^{†} c_{k} ρ_{b}^{e q}) = f (ϵ_{k})$ [where $f (ϵ_{k}) = \frac{1}{1 + e^{β ϵ_{k}}}$ is the fermi function], $β = \frac{1}{k_{B} T}$ ⁠, and

ρ_{s} (0) = ρ_{0} |ψ_{0}^{e l}⟩ ⟨ψ_{0}^{e l}| + ρ_{1} |ψ_{1}^{e l}⟩ ⟨ψ_{1}^{e l}| .

(7)

Note that ρ₀ + ρ₁ corresponds to one particular electronic state of the molecule so that ρ₀ and ρ₁ are $\in {H_{v i b} \otimes H_{d i f f} \otimes H_{s o l v e n t}}$ ⁠. Term No. 1b corresponding to H_c(t)H_c(t − τ)ρ(t) becomes

\begin{align} T r_{b} \sum_{k, k^{'}} (e^{- \frac{i ϵ_{k} t}{ℏ}} e^{\frac{i ϵ_{k^{'}} (t - τ)}{ℏ}} c_{k} c_{k^{'}}^{†} e^{\frac{i H_{s} t}{ℏ}} W_{k} (\vec{r}) d^{†} e^{- \frac{i H_{s} t}{ℏ}}) \\ (\times e^{\frac{i H_{s} (t - τ)}{ℏ}} W_{k^{'}} (\vec{r}) d e^{- \frac{i H_{s} (t - τ)}{ℏ}}) ρ_{b}^{e q} \otimes ρ_{s} (t) \\ = \sum_{k} e^{- \frac{i ϵ_{k} τ}{ℏ}} (1 - f (ϵ_{k})) e^{\frac{i H_{s} t}{ℏ}} W_{k} (\vec{r}) d e^{- \frac{i H_{s} τ}{ℏ}} W_{k} (\vec{r}) d^{†} e^{- \frac{i H_{s} (t - τ)}{ℏ}} \\ \times e^{\frac{i H_{s} t}{ℏ}} ρ_{s} (0) e^{- \frac{i H_{s} t}{ℏ}} \\ = \sum_{k} e^{- \frac{i ϵ_{k} τ}{ℏ}} (1 - f (ϵ_{k})) e^{\frac{i H_{s} t}{ℏ}} e^{- \frac{i H_{0} τ}{ℏ}} W_{k} (\vec{r}) e^{\frac{i H_{1} τ}{ℏ}} W_{k} (\vec{r}) ρ_{1} \\ \times |ψ_{1}^{e l}⟩ ⟨ψ_{1}^{e l}| e^{- \frac{i H_{s} t}{ℏ}} . \end{align}

(8)

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:

\begin{align} T r_{b} \sum_{k, k^{'}} e^{\frac{i ϵ_{k} t}{ℏ}} e^{- \frac{i ϵ_{k^{'}} (t - τ)}{ℏ}} (c_{k}^{†} e^{\frac{i H_{s} t}{ℏ}} W_{k} (\vec{r}) d e^{- \frac{i H_{s} t}{ℏ}}) (ρ_{b}^{e q} \otimes ρ_{s} (t)) \\ \times (c_{k^{'}} e^{\frac{i H_{s} (t - τ)}{ℏ}} W_{k^{'}} (\vec{r}) d^{†} e^{- \frac{i H_{s} (t - τ)}{ℏ}}) \\ = \sum_{k} e^{\frac{i ϵ_{k} τ}{ℏ}} (1 - f (ϵ_{k})) e^{\frac{i H_{s} t}{ℏ}} W_{k} (\vec{r}) d ρ_{s} (0) e^{- \frac{i H_{s} τ}{ℏ}} \\ \times W_{k} (\vec{r}) d^{†} e^{- \frac{i H_{s} (t - τ)}{ℏ}} \\ = \sum_{k} e^{\frac{i ϵ_{k} τ}{ℏ}} (1 - f (ϵ_{k})) e^{\frac{i H_{s} t}{ℏ}} W_{k} (\vec{r}) ρ_{1} |ψ_{0}^{e l}⟩ ⟨ψ_{0}^{e l}| e^{- \frac{i H_{1} τ}{ℏ}} \\ \times W_{k} (\vec{r}) e^{\frac{i H_{0} τ}{ℏ}} e^{- \frac{i H_{s} t}{ℏ}} \end{align}

(9)

and term No. 2b:

\begin{align} T r_{b} \sum_{k, k^{'}} e^{- \frac{i ϵ_{k} t}{ℏ}} e^{\frac{i ϵ_{k^{'}} (t - τ)}{ℏ}} (c_{k} e^{\frac{i H_{s} t}{ℏ}} W_{k} (\vec{r}) d^{†} e^{- \frac{i H_{s} t}{ℏ}}) (ρ_{b}^{e q} \otimes ρ_{s} (t)) \\ \times (c_{k^{'}}^{†} e^{\frac{i H_{s} (t - τ)}{ℏ}} W_{k^{'}} (\vec{r}) d e^{- \frac{i H_{s} (t - τ)}{ℏ}}) \\ = \sum_{k} e^{- \frac{i ϵ_{k} τ}{ℏ}} f (ϵ_{k}) e^{\frac{i H_{s} t}{ℏ}} W_{k} (\vec{r}) d^{†} ρ_{s} (0) e^{- \frac{i H_{s} τ}{ℏ}} W_{k} (\vec{r}) d e^{- \frac{i H_{s} (t - τ)}{ℏ}} \\ = \sum_{k} e^{- \frac{i ϵ_{k} τ}{ℏ}} f (ϵ_{k}) e^{\frac{i H_{s} t}{ℏ}} W_{k} (\vec{r}) ρ_{0} |ψ_{1}^{e l}⟩ ⟨ψ_{1}^{e l}| e^{- \frac{i H_{0} τ}{ℏ}} \\ \times W_{k} (\vec{r}) e^{\frac{i H_{1} τ}{ℏ}} e^{- \frac{i H_{s} t}{ℏ}} . \end{align}

(10)

The third set of terms, corresponding to ρ(t)H_c(t − τ)H_c(t), becomes

term No. 3a:

\begin{align} T r_{b} \sum_{k, k^{'}} e^{- \frac{i ϵ_{k} (t - τ)}{ℏ}} e^{\frac{i ϵ_{k^{'}} t}{ℏ}} (ρ_{b}^{e q} \otimes ρ_{s} (t)) (c_{k} e^{\frac{i H_{s} (t - τ)}{ℏ}} W_{k} (\vec{r}) d^{†} e^{- \frac{i H_{s} (t - τ)}{ℏ}}) \\ \times (c_{k^{'}}^{†} e^{\frac{i H_{s} t}{ℏ}} W_{k^{'}} (\vec{r}) d e^{- \frac{i H_{s} t}{ℏ}}) \\ = \sum_{k} e^{\frac{i ϵ_{k} τ}{ℏ}} (1 - f (ϵ_{k})) e^{\frac{i H_{s} t}{ℏ}} ρ_{s} (0) e^{- \frac{i H_{s} τ}{ℏ}} W_{k} (\vec{r}) d^{†} e^{\frac{i H_{s} τ}{ℏ}} \\ \times W_{k} (\vec{r}) d e^{- \frac{i H_{s} t}{ℏ}} \\ = \sum_{k} e^{\frac{i ϵ_{k} τ}{ℏ}} (1 - f (ϵ_{k})) e^{\frac{i H_{s} t}{ℏ}} ρ_{1} |ψ_{1}^{e l}⟩ ⟨ψ_{1}^{e l}| e^{- \frac{i H_{1} τ}{ℏ}} W_{k} (\vec{r}) e^{\frac{i H_{0} τ}{ℏ}} \\ \times W_{k} (\vec{r}) e^{- \frac{i H_{s} t}{ℏ}} \end{align}

(11)

and term No. 3b:

\begin{align} T r_{b} \sum_{k, k^{'}} e^{\frac{i ϵ_{k} (t - τ)}{ℏ}} e^{- \frac{i ϵ_{k^{'}} t}{ℏ}} (ρ_{b}^{e q} \otimes ρ_{s} (t)) (c_{k}^{†} e^{\frac{i H_{s} (t - τ)}{ℏ}} W_{k} (\vec{r}) d e^{- \frac{i H_{s} (t - τ)}{ℏ}}) \\ \times (c_{k^{'}} e^{\frac{i H_{s} t}{ℏ}} W_{k^{'}} (\vec{r}) d^{†} e^{- \frac{i H_{s} t}{ℏ}}) \\ = \sum_{k} e^{- \frac{i ϵ_{k} τ}{ℏ}} f (ϵ_{k}) e^{\frac{i H_{s} t}{ℏ}} ρ_{s} (0) e^{- \frac{i H_{s} τ}{ℏ}} W_{k} (\vec{r}) d e^{\frac{i H_{s} τ}{ℏ}} W_{k} (\vec{r}) d^{†} e^{- \frac{i H_{s} t}{ℏ}} \\ = \sum_{k} e^{- \frac{i ϵ_{k} τ}{ℏ}} f (ϵ_{k}) e^{\frac{i H_{s} t}{ℏ}} ρ_{0} |ψ_{0}^{e l}⟩ ⟨ψ_{0}^{e l}| e^{- \frac{i H_{0} τ}{ℏ}} W_{k} (\vec{r}) e^{\frac{i H_{1} τ}{ℏ}} \\ \times W_{k} (\vec{r}) e^{- \frac{i H_{s} t}{ℏ}} . \end{align}

(12)

Finally, the fourth set of term H_c(t − τ)ρ(t)H_c(t) yields

term No. 4a:

\begin{align} T r_{b} \sum_{k, k^{'}} e^{\frac{i ϵ_{k} (t - τ)}{ℏ}} e^{- \frac{i ϵ_{k^{'}} t}{ℏ}} (c_{k}^{†} e^{\frac{i H_{s} (t - τ)}{ℏ}} W_{k} (\vec{r}) d e^{- \frac{i H_{s} (t - τ)}{ℏ}}) (ρ_{b}^{e q} \otimes ρ_{s} (t)) \\ \times (c_{k^{'}} e^{\frac{i H_{s} t}{ℏ}} W_{k^{'}} (\vec{r}) d^{†} e^{- \frac{i H_{s} t}{ℏ}}) \\ = \sum_{k} e^{- \frac{i ϵ_{k} τ}{ℏ}} (1 - f (ϵ_{k})) e^{\frac{i H_{s} (t - τ)}{ℏ}} W_{k} (\vec{r}) d e^{\frac{i H_{s} τ}{ℏ}} ρ_{s} (0) \\ \times W_{k} (\vec{r}) d^{†} e^{- \frac{i H_{s} t}{ℏ}} \\ = \sum_{k} e^{- \frac{i ϵ_{k} τ}{ℏ}} (1 - f (ϵ_{k})) e^{\frac{i H_{s} t}{ℏ}} e^{- \frac{i H_{0} τ}{ℏ}} W_{k} (\vec{r}) e^{\frac{i H_{1} τ}{ℏ}} ρ_{1} |ψ_{0}^{e l}⟩ ⟨ψ_{0}^{e l}| \\ \times W_{k} (\vec{r}) e^{- \frac{i H_{s} t}{ℏ}} \end{align}

(13)

term No. 4b:

\begin{align} T r_{b} \sum_{k, k^{'}} e^{- \frac{i ϵ_{k} (t - τ)}{ℏ}} e^{\frac{i ϵ_{k^{'}} t}{ℏ}} (c_{k} e^{\frac{i H_{s} (t - τ)}{ℏ}} W_{k} (\vec{r}) d^{†} e^{- \frac{i H_{s} (t - τ)}{ℏ}}) (ρ_{b}^{e q} \otimes ρ_{s} (t)) \\ \times (c_{k^{'}}^{†} e^{\frac{i H_{s} t}{ℏ}} W_{k^{'}} (\vec{r}) d e^{- \frac{i H_{s} t}{ℏ}}) \\ = \sum_{k} e^{\frac{i ϵ_{k} τ}{ℏ}} f (ϵ_{k}) e^{\frac{i H_{s} (t - τ)}{ℏ}} W_{k} (\vec{r}) d^{†} e^{\frac{i H_{s} τ}{ℏ}} ρ_{s} (0) W_{k} (\vec{r}) d e^{- \frac{i H_{s} t}{ℏ}} \\ = \sum_{k} e^{\frac{i ϵ_{k} τ}{ℏ}} f (ϵ_{k}) e^{\frac{i H_{s} t}{ℏ}} e^{- \frac{i H_{1} τ}{ℏ}} W_{k} (\vec{r}) e^{\frac{i H_{0} τ}{ℏ}} ρ_{0} |ψ_{1}^{e l}⟩ ⟨ψ_{1}^{e l}| \\ \times W_{k} (\vec{r}) e^{- \frac{i H_{s} t}{ℏ}} . \end{align}

(14)

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 ρ₀ and ρ₁,

\begin{align} \frac{\partial ρ_{0}}{\partial t} & = - \frac{i}{ℏ} [H_{0}, ρ_{0}] - \sum_{k} \frac{1}{ℏ^{2}} \int_{0}^{\infty} d τ \\ \times (e^{- \frac{i ϵ_{k} τ}{ℏ}} f (ϵ_{k}) ρ_{0} e^{- \frac{i H_{0} τ}{ℏ}} W_{k} (\vec{r}) e^{\frac{i H_{1} τ}{ℏ}} W_{k} (\vec{r})) \\ - e^{- \frac{i ϵ_{k} τ}{ℏ}} (1 - f (ϵ_{k})) e^{- \frac{i H_{0} τ}{ℏ}} W_{k} (\vec{r}) e^{\frac{i H_{1} τ}{ℏ}} ρ_{1} W_{k} (\vec{r}) \\ - e^{\frac{i ϵ_{k} τ}{ℏ}} (1 - f (ϵ_{k})) W_{k} (\vec{r}) ρ_{1} e^{- \frac{i H_{1} τ}{ℏ}} W_{k} (\vec{r}) e^{\frac{i H_{0} τ}{ℏ}} \\ + (e^{\frac{i ϵ_{k} τ}{ℏ}} f (ϵ_{k}) W_{k} (\vec{r}) e^{- \frac{i H_{1} τ}{ℏ}} W_{k} (\vec{r}) e^{\frac{i H_{0} τ}{ℏ}} ρ_{0}) \end{align}

(15)

and

\begin{align} \frac{\partial ρ_{1}}{\partial t} & = - \frac{i}{ℏ} [H_{1}, ρ_{1}] - \sum_{k} \frac{1}{ℏ^{2}} \int_{0}^{\infty} d τ \\ \times (e^{\frac{i ϵ_{k} τ}{ℏ}} (1 - f (ϵ_{k})) ρ_{1} e^{- \frac{i H_{1} τ}{ℏ}} W_{k} (\vec{r}) e^{\frac{i H_{0} τ}{ℏ}} W_{k} (\vec{r})) \\ - e^{\frac{i ϵ_{k} τ}{ℏ}} f (ϵ_{k}) e^{- \frac{i H_{1} τ}{ℏ}} W_{k} (\vec{r}) e^{\frac{i H_{0} τ}{ℏ}} ρ_{0} W_{k} (\vec{r}) \\ - e^{- \frac{i ϵ_{k} τ}{ℏ}} f (ϵ_{k}) W_{k} (\vec{r}) ρ_{0} e^{- \frac{i H_{0} τ}{ℏ}} W_{k} (\vec{r}) e^{\frac{i H_{1} τ}{ℏ}} \\ + (e^{- \frac{i ϵ_{k} τ}{ℏ}} (1 - f (ϵ_{k})) W_{k} (\vec{r}) e^{- \frac{i H_{0} τ}{ℏ}} W_{k} (\vec{r}) e^{\frac{i H_{1} τ}{ℏ}} ρ_{1}) . \end{align}

(16)

These equations can be represented in the following form:

\frac{\partial ρ}{\partial t} = - \frac{i}{ℏ} [H, ρ] - \hat{\hat{L}} ρ,

(17)

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

\begin{aligned} \frac{\partial ρ_{0}}{\partial t} & = - \frac{i}{ℏ} [H_{0}, ρ_{0}] - \hat{\hat{L}} ρ_{0}, \\ \frac{\partial ρ_{1}}{\partial t} & = - \frac{i}{ℏ} [H_{1}, ρ_{1}] - \hat{\hat{L}} ρ_{1}, \end{aligned}

(18)

with ρ_s partitioned into ρ₀ and ρ₁ using Eq. (7). The superoperator $\hat{\hat{L}}$ is defined by its operation on ρ₀ and ρ₁,

\begin{align} \hat{\hat{L}} ρ_{0} & = \sum_{k} \frac{1}{ℏ^{2}} \int_{0}^{\infty} d τ e^{\frac{i ϵ_{k} τ}{ℏ}} f (ϵ_{k}) W_{k} (\vec{r}) e^{- \frac{i H_{1} τ}{ℏ}} W_{k} (\vec{r}) e^{\frac{i H_{0} τ}{ℏ}} ρ_{0} \\ - e^{\frac{i ϵ_{k} τ}{ℏ}} (1 - f (ϵ_{k})) W_{k} (\vec{r}) ρ_{1} e^{- \frac{i H_{1} τ}{ℏ}} W_{k} (\vec{r}) e^{\frac{i H_{0} τ}{ℏ}} + h . c . \end{align}

(19)

and

\begin{align} \hat{\hat{L}} ρ_{1} & = \sum_{k} \frac{1}{ℏ^{2}} \int_{0}^{\infty} d τ e^{- \frac{i ϵ_{k} τ}{ℏ}} (1 - f (ϵ_{k})) W_{k} (\vec{r}) e^{- \frac{i H_{0} τ}{ℏ}} W_{k} (\vec{r}) e^{\frac{i H_{1} τ}{ℏ}} ρ_{1} \\ - e^{- \frac{i ϵ_{k} τ}{ℏ}} f (ϵ_{k}) W_{k} (\vec{r}) ρ_{0} e^{- \frac{i H_{0} τ}{ℏ}} W_{k} (\vec{r}) e^{\frac{i H_{1} τ}{ℏ}} + h . c ., \end{align}

(20)

where both ρ₀ and ρ₁ 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 $ρ_{0}^{v i b} (x, ζ)$ and $ρ_{1}^{v i b} (x, ζ)$ ⁠. Note that $ρ_{0}^{v i b}$ and $ρ_{1}^{v i b}$ are $\in {H_{p r o t o n}}$ ⁠, 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 ${(A B)}_{W} = A_{W} e^{\frac{ℏ Λ}{2 i}} B_{W}$ ⁠, where Λ is the Poisson bracket operator $\overset{\leftrightarrow}{Λ} = {\overset{⃖}{\nabla}}_{P} \cdot {\vec{\nabla}}_{R} - {\overset{⃖}{\nabla}}_{R} \cdot {\vec{\nabla}}_{P}$ ⁠. Then, following the standard notation of classical mechanics, let us define ${Ô_{1}, Ô_{2}} = \sum_{α} \frac{\partial Ô_{1}}{\partial R_{α}} \frac{\partial Ô_{2}}{\partial p_{α}} - \frac{\partial Ô_{1}}{\partial p_{α}} \frac{\partial Ô_{2}}{\partial R_{α}}$ as the usual Poisson bracket. Note that, because Ô₁ and Ô₂ are operators, the order of the operators is critically important. Thus, $\frac{\partial Ô_{1}}{\partial x_{α}} \frac{\partial Ô_{2}}{\partial p_{α}} \neq \frac{\partial Ô_{1}}{\partial p_{α}} \frac{\partial Ô_{2}}{\partial x_{α}}$ ⁠. Finally, when we calculate the Wigner transformation of the commutator $- \frac{i}{ℏ} [H_{s}, ρ_{s}]$ to first order in ℏ, which is a classical assumption about the diffusion (x) and solvent (ζ) coordinates (Assumption No. 4), we find the result $- \frac{i}{ℏ} [H_{s}, ρ_{s}] = - \frac{i}{ℏ} [H_{s}^{v i b}, ρ_{s}^{v i b}] + {H_{s}^{v i b}, ρ_{s}^{v i b}}_{a}$ ⁠, where ${Ô_{1}, Ô_{2}}_{a} = \frac{1}{2} ({Ô_{1}, Ô_{2}} - {Ô_{2}, Ô_{1}})$ ⁠. We are left with the famous quantum classical Liouville equation (QCLE) for ${\hat{ρ}}_{s}^{v i b}$ ⁠,

\begin{align} \begin{aligned} \frac{\partial {\hat{ρ}}_{0}^{v i b}}{\partial t} & = \frac{1}{2} ({Ĥ_{0}^{v i b}, {\hat{ρ}}_{0}^{v i b}} - {{\hat{ρ}}_{0}^{v i b}, Ĥ_{0}^{v i b}}) - \frac{i}{ℏ} [Ĥ_{0}^{v i b}, {\hat{ρ}}_{0}^{v i b}] - {\hat{\hat{L}}}^{v i b} {\hat{ρ}}_{0}^{v i b}, \\ \frac{\partial {\hat{ρ}}_{1}^{v i b}}{\partial t} & = \frac{1}{2} ({Ĥ_{1}^{v i b}, {\hat{ρ}}_{1}^{v i b}} - {{\hat{ρ}}_{1}^{v i b}, Ĥ_{1}^{v i b}}) - \frac{i}{ℏ} [Ĥ_{1}^{v i b}, {\hat{ρ}}_{1}^{v i b}] - {\hat{\hat{L}}}^{v i b} {\hat{ρ}}_{1}^{v i b} . \end{aligned} \end{align}

(21)

One would like to assume that, just as for the QCLE,³⁰ Eq. (21) is correct to order $O (ℏ)$ ⁠. However, in Eq. (21), we have implicitly made one more (fairly standard) assumption, i.e., that ${(\hat{\hat{L}} ρ_{s})}^{v i b} = {\hat{\hat{L}}}^{v i b} {\hat{ρ_{s}}}^{v i b}$ ⁠,³¹ 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 ${(\hat{\hat{L}} ρ_{s})}^{v i b}$ ⁠. 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 ${\hat{\hat{L}}}^{v i b} {\hat{ρ}}_{0}^{v i b}$ in vibrational states μ, ν [using Eqs. (19) and (20)],

\begin{align} ⟨μ| {\hat{\hat{L}}}^{v i b} {\hat{ρ}}_{0}^{v i b} |ν⟩ & = - \sum_{k, μ^{'}, ν^{'}, ξ} \frac{{|W_{k}|}^{2}}{ℏ^{2}} \int_{0}^{\infty} d τ (e^{- \frac{i ϵ_{k} τ}{ℏ}} f (ϵ_{k}) ⟨μ| ρ_{0} |ξ⟩) \\ (\times ⟨ξ| e^{- \frac{i H_{0} τ}{ℏ}} |μ^{'}⟩ ⟨μ^{'}| e^{\frac{i H_{1} τ}{ℏ}} |ν^{'}⟩) \\ - e^{- \frac{i ϵ_{k} τ}{ℏ}} (1 - f (ϵ_{k})) ⟨μ| e^{- \frac{i H_{0} τ}{ℏ}} |μ^{'}⟩ ⟨μ^{'}| e^{\frac{i H_{1} τ}{ℏ}} ρ_{1} |ν^{'}⟩ ⟨ν^{'}| |ν⟩ \\ - e^{\frac{i ϵ_{k} τ}{ℏ}} (1 - f (ϵ_{k})) ⟨μ| |μ^{'}⟩ ⟨μ^{'}| ρ_{1} |ν^{'}⟩ ⟨ν^{'}| e^{- \frac{i H_{1} τ}{ℏ}} e^{\frac{i H_{0} τ}{ℏ}} |ν⟩ \\ + (e^{\frac{i ϵ_{k} τ}{ℏ}} f (ϵ_{k}) ⟨μ| |μ^{'}⟩ ⟨μ^{'}| e^{- \frac{i H_{1} τ}{ℏ}} |ξ⟩ ⟨ξ| e^{\frac{i H_{0} τ}{ℏ}} ρ_{0} |ν⟩), \end{align}

(22)

where primed (unprimed) Greek indices correspond to proton vibrational states in the basis of electronic state 1 (state 0). We can similarly expand ${\hat{\hat{L}}}^{v i b} {\hat{ρ}}_{1}^{v i b}$ in vibrational states μ′, ν′ to obtain

\begin{align} ⟨μ^{'}| {\hat{\hat{L}}}^{v i b} {\hat{ρ}}_{1}^{v i b} |ν^{'}⟩ & = - \sum_{k, μ, ν, ξ^{'}} \frac{{|W_{k}|}^{2}}{ℏ^{2}} \int_{0}^{\infty} d τ (e^{\frac{i ϵ_{k} τ}{ℏ}} (1 - f (ϵ_{k})) ⟨μ^{'}| ρ_{1} |ξ^{'}⟩) \\ (\times ⟨ξ^{'}| e^{- \frac{i H_{1} τ}{ℏ}} |μ⟩ ⟨μ| e^{\frac{i H_{0} τ}{ℏ}} |ν⟩) \\ - e^{\frac{i ϵ_{k} τ}{ℏ}} f (ϵ_{k}) ⟨μ^{'}| e^{- \frac{i H_{1} τ}{ℏ}} |μ⟩ ⟨μ| e^{\frac{i H_{0} τ}{ℏ}} ρ_{0} |ν⟩ ⟨ν| |ν^{'}⟩ \\ - e^{- \frac{i ϵ_{k} τ}{ℏ}} f (ϵ_{k}) ⟨μ^{'}| |μ⟩ ⟨μ| ρ_{0} |ν⟩ ⟨ν| e^{- \frac{i H_{0} τ}{ℏ}} e^{\frac{i H_{1} τ}{ℏ}} |ν^{'}⟩ \\ + (e^{- \frac{i ϵ_{k} τ}{ℏ}} (1 - f (ϵ_{k})) ⟨μ^{'}| |μ⟩ ⟨μ| e^{- \frac{i H_{0} τ}{ℏ}} |ξ^{'}⟩) \\ (\times ⟨ξ^{'}| e^{\frac{i H_{1} τ}{ℏ}} ρ_{1} |ν^{'}⟩) . \end{align}

(23)

At this point, if we make the substitution $\int_{0}^{\infty} d τ e^{\frac{i ϵ_{k} τ}{ℏ}} = π δ (ϵ_{k}) ℏ$ ⁠, we arrive at a relatively complete expression for $\frac{\partial ρ_{0}^{μ ν}}{\partial t}$ ⁠,

\begin{align} \frac{\partial ρ_{0}^{μ ν}}{\partial t} & = \frac{1}{2} ⟨μ| ({Ĥ_{0}^{v i b}, {\hat{ρ}}_{0}^{v i b}} - {{\hat{ρ}}_{0}^{v i b}, Ĥ_{0}^{v i b}}) |ν⟩ - \frac{i}{ℏ} (U_{0}^{μ} - U_{0}^{ν}) ρ_{0}^{μ ν} \\ - \frac{Γ}{2 ℏ} \sum_{μ^{'}, ξ} f (U_{1}^{μ^{'}} - U_{0}^{ξ}) ⟨μ| |μ^{'}⟩ ⟨μ^{'}| |ξ⟩ ρ_{0}^{ξ ν} \\ - \frac{Γ}{2 ℏ} \sum_{μ^{'}, ξ} ρ_{0}^{μ ξ} f (U_{1}^{μ^{'}} - U_{0}^{ξ}) ⟨ξ| |μ^{'}⟩ ⟨μ^{'}| |ν⟩ \\ + \frac{Γ}{2 ℏ} \sum_{μ^{'}, ν^{'}} (1 - f (U_{1}^{ν^{'}} - U_{0}^{ν})) ⟨μ| |μ^{'}⟩ ⟨ν^{'}| |ν⟩ ρ_{1}^{μ^{'} ν^{'}} \\ + \frac{Γ}{2 ℏ} \sum_{μ^{'}, ν^{'}} ρ_{1}^{μ^{'} ν^{'}} (1 - f (U_{1}^{μ^{'}} - U_{0}^{μ})) ⟨μ| |μ^{'}⟩ ⟨ν^{'}| |ν⟩, \end{align}

(24)

and similarly for $\frac{\partial ρ_{1}^{μ^{'} ν^{'}}}{\partial t}$ ⁠,

\begin{align} \frac{\partial ρ_{1}^{μ^{'} ν^{'}}}{\partial t} & = \frac{1}{2} ⟨μ^{'}| ({Ĥ_{1}^{v i b}, {\hat{ρ}}_{1}^{v i b}} - {{\hat{ρ}}_{1}^{v i b}, Ĥ_{1}^{v i b}}) |ν^{'}⟩ - \frac{i}{ℏ} (U_{1}^{μ^{'}} - U_{1}^{ν^{'}}) ρ_{1}^{μ^{'} ν^{'}} \\ - \frac{Γ}{2 ℏ} \sum_{μ, ξ^{'}} (1 - f (U_{1}^{μ^{'}} - U_{0}^{ξ})) ⟨μ^{'}| |μ⟩ ⟨μ| |ξ^{'}⟩ ρ_{1}^{ξ^{'} ν^{'}} \\ - \frac{Γ}{2 ℏ} \sum_{μ, ξ^{'}} ρ_{1}^{μ^{'} ξ^{'}} (1 - f (U_{1}^{μ^{'}} - U_{0}^{ξ})) ⟨ξ^{'}| |μ⟩ ⟨μ| |ν^{'}⟩ \\ + \frac{Γ}{2 ℏ} \sum_{μ, ν} f (U_{1}^{ν^{'}} - U_{0}^{ν}) ⟨μ^{'}| |μ⟩ ⟨ν| |ν^{'}⟩ ρ_{0}^{μ ν} \\ + \frac{Γ}{2 ℏ} \sum_{μ, ν} ρ_{0}^{μ ν} f (U_{1}^{μ^{'}} - U_{0}^{μ}) ⟨μ^{'}| |μ⟩ ⟨ν| |ν^{'}⟩ . \end{align}

(25)

Here, $ρ_{0}^{μ ν} = ⟨μ| {\hat{ρ}}_{0}^{v i b} |ν⟩$ and $Γ = 2 π \sum_{k} {|W_{k}|}^{2} δ (ϵ - ϵ_{k})$ is the metal–molecule coupling in the wide-band limit that characterizes the rate with which electrons can hop between molecule and metal.

Equations (24) and (25) capture the coupling of electron transfer to classical solvent motion as well as to quantum mechanical proton motion. In principle, these equations can be solved using an advanced surface hopping algorithm with an unconstrained number of electrons.²⁹

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

\begin{aligned} \frac{\partial ρ_{0}^{μ μ}}{\partial t} & = ⟨μ| ({Ĥ_{0}^{v i b}, {\hat{ρ}}_{0}^{v i b}}) |μ⟩ - \frac{Γ}{ℏ} \sum_{μ^{'}} {|F_{μ \to μ^{'}}|}^{2} f (U_{1}^{μ^{'}} - U_{0}^{μ}) ρ_{0}^{μ μ} \\ + \frac{Γ}{ℏ} \sum_{μ^{'}} {|F_{μ \to μ^{'}}|}^{2} (1 - f (U_{1}^{μ^{'}} - U_{0}^{μ})) ρ_{1}^{μ^{'} μ^{'}} \\ \frac{\partial ρ_{1}^{μ^{'} μ^{'}}}{\partial t} & = ⟨μ^{'}| ({Ĥ_{1}^{v i b}, {\hat{ρ}}_{1}^{v i b}}) |μ^{'}⟩ + \frac{Γ}{ℏ} \sum_{μ} {|F_{μ \to μ^{'}}|}^{2} f (U_{1}^{μ^{'}} - U_{0}^{μ}) ρ_{0}^{μ μ} \\ - \frac{Γ}{ℏ} \sum_{μ} {|F_{μ \to μ^{'}}|}^{2} (1 - f (U_{1}^{μ^{'}} - U_{0}^{μ})) ρ_{1}^{μ^{'} μ^{'}}, \end{aligned}

(26)

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.³⁵

6. Final equations

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

\begin{aligned} \frac{\partial ρ_{0}^{μ}}{\partial t} & = {({H_{0}, ρ_{0}})}^{μ} - \sum_{μ^{'}} (γ_{0 \to 1}^{μ μ^{'}} ρ_{0}^{μ} - γ_{1 \to 0}^{μ μ^{'}} ρ_{1}^{μ^{'}}), \\ \frac{\partial ρ_{1}^{μ^{'}}}{\partial t} & = {({H_{1}, ρ_{1}})}^{μ^{'}} - \sum_{μ} (γ_{1 \to 0}^{μ μ^{'}} ρ_{1}^{μ^{'}} - γ_{0 \to 1}^{μ μ^{'}} ρ_{0}^{μ}), \end{aligned}

(27)

where the hopping rates are

\begin{aligned} γ_{0 \to 1}^{μ μ^{'}} & = \frac{Γ}{ℏ} f (U_{1}^{μ^{'}} - U_{0}^{μ}) {|F_{μ \to μ^{'}}|}^{2}, \\ γ_{1 \to 0}^{μ μ^{'}} & = \frac{Γ}{ℏ} (1 - f (U_{1}^{μ^{'}} - U_{0}^{μ})) {|F_{μ \to μ^{'}}|}^{2} . \end{aligned}

(28)

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, $U_{1}^{μ^{'}} - U_{0}^{μ}$ ⁠, instead of the more straightforward electronic diabatic energy difference U₁ − U₀. 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 $U_{0}^{μ} (x, ζ)$ and $U_{1}^{μ} (x, ζ)$ ⁠, which will be recast as $U_{A}^{μ}$ and $U_{B}^{μ^{'}}$ (μ and μ′ index proton vibrational levels). This naming convention is to allow for more natural comparison to the chemical reaction occurring at the electrode,

A + e^{-} \overset{k_{f}}{\underset{k_{b}}{⇌}} B,

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,

U_{A}^{μ} (ζ) = \frac{1}{2} m ω^{2} {(ζ_{A} - ζ)}^{2} + Ũ_{A}^{μ},

(29a)

U_{B}^{μ^{'}} (ζ) = \frac{1}{2} m ω^{2} {(ζ_{B} - ζ)}^{2} + Δ G + Ũ_{B}^{μ^{'}},

(29b)

\begin{align} E^{μ, μ^{'}} (ζ) & = U_{B}^{μ^{'}} - U_{A}^{μ} \\ = \frac{1}{2} m ω^{2} (ζ_{B}^{2} - ζ_{A}^{2} - 2 (ζ_{B} - ζ_{A}) ζ) + Δ G + (Ũ_{B}^{μ^{'}} - Ũ_{A}^{μ}) . \end{align}

(29c)

Here, $U_{A}^{μ}$ and $U_{B}^{μ^{'}}$ correspond to the vibronic contribution to the energy for electronic states A and B and proton vibrational states μ and μ′, respectively. $Ũ_{A}^{μ}$ and $Ũ_{B}^{μ^{'}}$ 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 $U_{A}^{μ} (ζ)$ and $U_{B}^{μ} (ζ)$ ⁠, respectively. Therefore, the solvent reorganization energy can be defined as $λ = \frac{1}{2} m ω^{2} {(ζ_{B} - ζ_{A})}^{2}$ ⁠.^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 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:

\frac{\partial c_{A}^{μ}}{\partial t} = D_{A}^{x} \frac{\partial^{2} c_{A}^{j}}{\partial x^{2}} - \sum_{μ^{'}} (k_{f}^{μ \to μ^{'}} c_{A}^{μ} - k_{b}^{μ^{'} \to μ} c_{B}^{μ^{'}}) δ (x),

(30a)

\frac{\partial c_{B}^{μ^{'}}}{\partial t} = D_{B}^{x} \frac{\partial^{2} c_{B}^{j}}{\partial x^{2}} + \sum_{μ} (k_{f}^{μ \to μ^{'}} c_{A}^{μ} - k_{b}^{μ^{'} \to μ} c_{B}^{μ^{'}}) δ (x),

(30b)

where μ, μ′ is an index for the proton vibrational state in species A and B, respectively, and $k_{f / b}^{μ \to μ^{'}}$ are the forward/backward hopping rates between vibrational states μ and μ′. $c_{A}^{μ}$ and $c_{B}^{μ^{'}}$ 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 (⁠ $D_{A}^{x} = D_{B}^{x} \equiv D$ ⁠). To visualize these combined electronic-vibrational energy surfaces, see Fig. 1.

FIG. 1.

View large Download slide

Diabatic PESs for A and B for a system with three proton vibrational states, where μ is the vibrational state index. Here, the proton vibrational states are treated as harmonic oscillator states, with energy spacing ℏω_p. All parameters are in atomic units (a.u.) and are m = 2000, ω = 0.000 25 [from Eq. (29)], ω_p = 0.012, Γ = 10⁻⁵, ΔG = 0, and kT = 0.000 95.

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 $k_{f}^{μ \to μ^{'}}$ and $k_{b}^{μ^{'} \to μ}$ [see Eq. (32) below]. For all results below, we will choose spatial and initial boundary conditions as follows:

c_{A}^{μ} (t = 0, x) = \frac{e^{- β U_{A}^{μ}}}{\sum_{ν} e^{- β U_{A}^{ν}}},

(31a)

c_{B}^{μ} (t = 0, x) = 0,

(31b)

c_{A}^{μ} (t, x = - 1) = c_{A}^{μ} (t, x = 0),

(31c)

c_{B}^{μ} (t, x = - 1) = c_{B}^{μ} (t, x = 0),

(31d)

c_{A}^{μ} (t, x = \infty) = \frac{e^{- β U_{A}^{μ}}}{\sum_{ν} e^{- β U_{A}^{ν}}},

(31e)

c_{B}^{μ} (t, x = \infty) = 0 .

(31f)

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.³⁷

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,

k_{f}^{μ \to μ^{'}} = \frac{Γ}{ℏ} \frac{\int_{- \infty}^{\infty} d ζ e^{- β U_{A}^{μ} (ζ)} f (E^{μ, μ^{'}} (ζ))}{\int_{- \infty}^{\infty} d ζ e^{- β U_{A}^{μ} (ζ)}},

(32a)

k_{b}^{μ^{'} \to μ} = \frac{Γ}{ℏ} \frac{\int_{- \infty}^{\infty} d ζ e^{- β U_{B}^{μ^{'}} (ζ)} (1 - f (E^{μ, μ^{'}} (ζ)))}{\int_{- \infty}^{\infty} d ζ e^{- β U_{B}^{μ^{'}} (ζ)}} .

(32b)

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 ν,

Δ G (t) = Δ G_{i} + ν t,

(33)

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,

I (t) = \sum_{μ, μ^{'}} (k_{f}^{μ \to μ^{'}} (t) c_{A}^{μ} (t) - k_{b}^{μ^{'} \to μ} (t) c_{B}^{μ^{'}} (t)) |_{x = 0},

(34)

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,

I (t) = D_{A}^{x} \sum_{μ} \frac{\partial c_{A}^{μ}}{\partial x} |_{x = 0} = - D_{B}^{x} \sum_{μ^{'}} \frac{\partial c_{B}^{μ^{'}}}{\partial x} |_{x = 0} .

(35)

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³¹ and, more recently, to study electrochemical ET under liner sweep conditions.²⁵ 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:

\frac{\partial c_{A}^{μ}}{\partial t} = L_{A} c_{A}^{μ} - \sum_{μ^{'}} (k_{f}^{μ \to μ^{'}} c_{A}^{μ} - k_{b}^{μ^{'} \to μ} c_{B}^{μ^{'}}) δ (x),

(36a)

\frac{\partial c_{B}^{μ^{'}}}{\partial t} = L_{B} c_{B}^{μ^{'}} + \sum_{μ} (k_{f}^{μ \to μ^{'}} c_{A}^{μ} - k_{b}^{μ^{'} \to μ} c_{B}^{μ^{'}}) δ (x) .

(36b)

Here, we define superoperators,

L_{A} = \frac{\partial}{\partial x} (D_{A} \frac{\partial}{\partial x}),

(37a)

L_{B} = \frac{\partial}{\partial x} (D_{B} \frac{\partial}{\partial x}) .

(37b)

To be as specific as possible, let us now consider the case of a two vibrational state system, with energy spacing ℏω_p. Combining $L_{A}, L_{B}, k_{f}$ ⁠, and k_b, we can obtain the general matrix equation for c_A and c_B,

\begin{aligned} (\begin{matrix} \dot{c_{A}^{0}} \\ \dot{c_{A}^{1}} \\ \dot{c_{B}^{0}} \\ \dot{c_{B}^{1}} \end{matrix}) & = (\begin{matrix} L_{A} - \sum_{i}^{2} k_{f}^{0 \to i} δ (x) & 0 & k_{b}^{0 \to 0} δ (x) & k_{b}^{1 \to 0} δ (x) \\ 0 & L_{A} - \sum_{i}^{2} k_{f}^{1 \to i} δ (x) & k_{b}^{0 \to 1} δ (x) & k_{b}^{1 \to 1} δ (x) \\ k_{f}^{0 \to 0} δ (x) & k_{f}^{1 \to 0} δ (x) & L_{B} - \sum_{i}^{2} k_{b}^{0 \to i} δ (x) & 0 \\ k_{f}^{0 \to 1} δ (x) & k_{f}^{1 \to 1} δ (x) & 0 & L_{B} - \sum_{i}^{2} k_{b}^{1 \to i} δ (x) \end{matrix}) (\begin{matrix} c_{A}^{0} \\ c_{A}^{1} \\ c_{B}^{0} \\ c_{B}^{1} \end{matrix}) + \vec{y} \\ \dot{\vec{c}} & = M \vec{c} + \vec{y} . \end{aligned}

(38)

The generalization to many vibrational states is straightforward. Here, in Eq. (38), the $\vec{y}$ 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., $\dot{\vec{c}} (x = \infty) = \vec{0}$ ⁠. In other words, $\vec{y} (x \neq \infty) = 0$ and $\vec{y} (x = \infty) = - (M \vec{c}) (x = \infty)$ ⁠. 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 $\vec{y} = 0 \to \vec{y} = - (M \vec{c})$ at x = ∞.

$M$ can be non-Hermitian, since the off-diagonal elements of $M$ are not equal (⁠ $k_{f}^{i \to i^{'}} \neq k_{b}^{i^{'} \to i}$ ⁠) 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

S = (\begin{matrix} 1 & 0 & 0 & 0 \\ 0 & e^{\frac{β ℏ ω_{p}}{2}} & 0 & 0 \\ 0 & 0 & e^{\frac{β Δ G}{2}} & 0 \\ 0 & 0 & 0 & e^{\frac{β (Δ G + ℏ ω_{p})}{2}} \end{matrix})

(39)

such that the new operator $\tilde{M}$ satisfies

\tilde{M} = S M S^{- 1} = (\begin{matrix} L_{A} - \sum_{i}^{2} k_{f}^{0 \to i} δ (x) & 0 & k_{b}^{0 \to 0} δ (x) e^{- \frac{β E^{0,0}}{2}} & k_{b}^{1 \to 0} δ (x) e^{- \frac{β E^{0,1}}{2}} \\ 0 & L_{A} - \sum_{i}^{2} k_{f}^{1 \to i} δ (x) & k_{b}^{0 \to 1} δ (x) e^{- \frac{β E^{1,0}}{2}} & k_{b}^{1 \to 1} δ (x) e^{- \frac{β E^{1,1}}{2}} \\ k_{f}^{0 \to 0} δ (x) e^{\frac{β E^{0,0}}{2}} & k_{f}^{1 \to 0} δ (x) e^{\frac{β E^{1,0}}{2}} & L_{B} - \sum_{i}^{2} k_{b}^{0 \to i} δ (x) & 0 \\ k_{f}^{0 \to 1} δ (x) e^{\frac{β E^{0,1}}{2}} & k_{f}^{1 \to 1} δ (x) e^{\frac{β E^{1,1}}{2}} & 0 & L_{B} - \sum_{i}^{2} k_{b}^{1 \to i} δ (x) . \end{matrix})

(40)

Note that $\tilde{M}$ is Hermitian and can be easily diagonalized (also that $S L_{A} S^{- 1} = L_{A}$ because D_A is a constant, with no x-dependence). While the eigenvectors of $\tilde{M}$ will differ from those of $M$ ⁠, $M$ and $\tilde{M}$ 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 $\tilde{M}$ at each time step, we can simulate a linear sweep voltammetry experiment. Once we have the eigenvectors P and eigenvalues Λ of $\tilde{M} = P Λ P^{T}$ ⁠, we can change variables and recast the differential equation for the time evolution of c_A and c_B into renormalized coordinates,

\vec{b} \equiv P^{T} S \vec{c},

(41a)

\vec{z} \equiv P^{T} S \vec{y} .

(41b)

The final equation of motion becomes

\dot{\vec{b}} = Λ \vec{b} + \vec{z} .

(41c)

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

\vec{b} (t) = {\vec{b}}_{e q} + (\vec{b} - {\vec{b}}_{e q}) e^{Λ t},

(42a)

{\vec{b}}_{e q} = - Λ^{- 1} \vec{z} .

(42b)

With Eq. (41a), we can change $\vec{b} (t)$ to $\vec{c} (t) = (\begin{matrix} c_{A} \\ c_{B} \end{matrix})$ using P and S: $\vec{c} = S^{- 1} P \vec{b}$ ⁠. 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⁻¹⁶ a.u.), Γ = 0.000 27 eV (10⁻⁵ a.u.), D = 5.5 · 10⁻⁴ cm²/s (4.75 · 10⁻⁴ a.u.), ω = 1.03 · 10¹³ s⁻¹ (0.000 25 a.u.), m = 1.82 · 10⁻²⁷ 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 $L_{x} = 6 \sqrt{D t_{T}}$ ⁠, where t_T is the total time of the simulation, defined by the starting and ending ΔG and ν, $t_{T} = \frac{| Δ G_{f} - Δ G_{i} |}{ν}$ ⁠, and the grid spacing is uniform such that $d x = \frac{L_{x}}{N_{x}}$ ⁠. 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 $Λ = \frac{k^{0}}{\sqrt{D ν \frac{F}{R T}}}$ ⁠.^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⁰. 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⁰, 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 $U_{0}^{μ}$ and $U_{1}^{μ^{'}}$ ⁠. 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₀ is the potential for the proton assuming the molecule is uncharged and V₁ is the potential for the proton assuming the molecule is charged, we can write

\begin{aligned} V_{0} (q) = \frac{1}{2} m ω_{p}^{2} q^{2}, \\ V_{1} (q) = \frac{1}{2} m ω_{p}^{2} {(q - g_{p})}^{2}, \end{aligned}

(43)

which allows for an analytical calculation of the Franck–Condon factors. Note that $[\frac{Π^{2}}{2 m_{q}} + V_{0} (q)] |μ⟩ = Ũ_{A}^{μ} |μ⟩$ ⁠. The Franck–Condon factors are

F_{μ \to μ^{'}} = ⟨μ| |μ^{'}⟩ = \int d q ϕ_{μ} (q + g_{p}) ϕ_{μ^{'}} (q),

(44)

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

F_{μ \to μ^{'}} = {\frac{p!}{Q!}}^{0.5} {(\frac{g_{p}}{\sqrt{2}})}^{Q - p} e^{- \frac{g_{p}^{2}}{4}} L_{p}^{Q - p} (\frac{g_{p}^{2}}{2}) {[s g n (μ^{'} - μ)]}^{μ - μ^{'}} .

(45)

Here, p/Q are the minimum/maximum of μ and μ′, sgn is the sign function, and $L_{m}^{n}$ is a generalized Laguerre polynomial. The reorganization energy in the proton coordinate, λ_p, is defined as $λ_{p} = \frac{ℏ ω_{p} g_{p}^{2}}{2}$ ⁠. 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 $\frac{I_{p}^{P C E T}}{I_{p}^{E T}}$ and (iii) the deviation of the peak potential between PCET and ET, $| E_{p}^{P C E T} - E_{p}^{E T} |$ ⁠. 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. $E_{p}^{E T}$ and $I_{p}^{E T}$ 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 $| E_{p}^{P C E T} - E_{p}^{E T} |$ 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⁰, through the incomplete orbital overlap of the two ground proton vibrational states in electronic states A and B [see Eq. (45)].

FIG. 2.

View large Download slide

[(a)–(d)] 1D Linear sweep voltammograms for the harmonic oscillator model at different values of ω_p and g_p. (a) ω_p = 8.42kT, (b) ω_p = 10.53kT, (c) ω_p = 12.63kT, and (d) ω_p = 14.74kT. (e) and (g) show the ratio of the peak current in the PCET simulation to that of the ET simulation, with respect to frequency ω_p and displacement in q (g_p), respectively. (f) and (h) show the deviation of the peak potential in the PCET simulation and the ET simulation with respect to frequency ω_p and g_p, respectively. In (g)–(h), the colors of points and lines correspond to their respective ω_p values from (a)–(d). The different curves correspond to different values of g_p ranging from 0.75 Å to 1.65 Å. The black trace shows the case with no proton coordinate for reference. From the data in (e)–(g), we see that there is no dependence on ω_p for either I_p or E_p, suggesting that all hops are between ground vibrational states in electronic states A and B. Thus, the dependence of I_p and E_p on g_p suggests that the main influence of PCET on the rate of ET is through the orbital overlaps between the two ground vibrational states at a given value of g_p, and therefore, the overall effect of the proton is just to reduce the value of the standard rate constant. For these systems, Λ values are well defined and range from Λ = 520.26 for g_p = 0.75 Å to Λ = 11.18 for g_p = 1.65 Å. Here, we are operating largely in the reversible regime.

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.⁴² In this instance, the PESs in the proton coordinate are

\begin{aligned} V_{0} (q) = D_{e} {(1 - e^{a q})}^{2}, \\ V_{1} (q) = D_{e} {(1 - e^{- a (q - g_{p})})}^{2}, \end{aligned}

(46)

where $a = \sqrt{\frac{m ω_{p}^{2}}{2 D_{e}}}$ ⁠, and we allow g_p to vary. The energy levels are given by $E_{n} = ℏ ω_{p} \frac{2 D_{e} + 4 n D_{e} - ℏ ω_{p} (n^{2} + n + \frac{1}{4})}{4 D_{e}}$ ⁠, 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.

FIG. 3.

View large Download slide

Example of V₀(q) and V₁(q) for the Morse potential, ω_p = 12.63kT, g_p = 0.741 Å.

FIG. 4.

View large Download slide

[(a) and (b)] 1D Linear sweep voltammograms for the Morse oscillator model at different values of ω_p and g_p. (a) ω_p = 10.53kT and (b) ω_p = 14.74kT. (c) Ratio of the peak current in the PCET simulation to that of the ET simulation with respect to displacement in q (g_p) of the two electronic states 0 and 1. (d) Deviation of the peak potential in the PCET simulation vs the ET simulation with respect to g_p. The black trace shows the case with no proton coordinate for reference. Qualitatively, we see some of the same behavior as the harmonic oscillator case, such as E_p increasing; although bizarrely, we find that I_p seems to modestly increase with g_p. Even more interestingly, at large g_p, we see that there is now a ω_p dependence on the peak position. This suggests that the mechanism of PCET is different from that of the harmonic oscillator proton PES, as hops can now occur between ground and excited states in 0 and 1, respectively, instead of between the ground state and ground state. Thus, Λ values for these reactions are not entirely well-defined. Using the 0–0 (ground-to-ground) transition as reference, the Λ values for these results range from Λ = 0.0541 for g_p = 0.63 Å to Λ = 4.38 · 10⁻⁷ for g_p = 0.95 Å. If the Λ values above are relevant, then, here, we operate in the quasi-reversible and irreversible regimes (whereas Fig. 2 operates in the reversible regime) (Fig. 3).

C. Quartic potential

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

V_{A} (q) = E_{r e f} (\frac{a_{4}^{A}}{L^{4}} {(q - q_{c})}^{4} - \frac{a_{2}^{A}}{L^{2}} {(q - q_{c})}^{2} + \frac{a_{1}^{A}}{L} (q - q_{c})),

(47a)

V_{B} (q) = E_{r e f} (\frac{a_{4}^{B}}{L^{4}} {(q_{c} - q)}^{4} - \frac{a_{2}^{B}}{L^{2}} {(q_{c} - q)}^{2} + \frac{a_{1}^{B}}{L} (q_{c} - q)),

(47b)

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., $a_{i}^{A} = a_{i}^{B}$ for i = 1, 2, 4). L is a reference length scale, $L = \sqrt{\frac{ℏ}{m ω_{p}}}$ ⁠, and E_ref is a reference energy scale. q_c is chosen as $\sqrt{\frac{a_{2} ℏ}{4 a_{4} m ω_{p}}}$ ⁠, and shown in Fig. 5 are example quartic PESs for four different values of a₁. When a₄ = a₁ = 0, a₂ = −1, this model is identical to the quadratic harmonic oscillator model. The parameter a₁ dictates the asymmetry of the two wells in the quartic model and yields symmetric double wells when a₁ = 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₁. When a₁ = 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₁ 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.⁴⁶ 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.

FIG. 5.

View large Download slide

Proton PESs for A and B in q for different values of a₁. m = 1.82 · 10⁻²⁷ kg, ω_p = 8.42kT, a₄ = 0.005, and a₂ = 0.05. Note that the PESs for A and B are symmetric with $q_{c} = \sqrt{\frac{2.5 ℏ}{m ω_{p}}}$ ⁠.

FIG. 6.

View large Download slide

1D Linear sweep voltammograms at different values of ω_p and a₁ for the quartic PES in the proton coordinate. (a) L = 0.25 (ω_p = 8.42kT) and E_ref = 8.42kT. (b) L = 0.204 (ω_p = 12.63kT) and E_ref = 12.63kT. a₄ = 0.005 and a₂ = 0.05. a₄ and a₂ are chosen such that q_c is the same for both models in (a) and (b). The black trace shows the case with no proton coordinate, for reference. (c) Ratio of the peak current in the PCET simulation to that of the ET simulation with respect to the parameter a₁ of the quartic PES. (d) Deviation of the peak potential in the PCET simulation and the ET simulation with respect to the parameter a₁ of the quartic PES. Qualitatively, we see very similar behavior to the anharmonic Morse oscillator. This finding reaffirms the notion that anharmonicity in the proton coordinate PES due to a₁ can give rise to slightly different effects in the voltammetric simulation, suggesting that knowledge of the proton PES is essential for modeling an IV experiment quantitatively accurately. As in Fig. 4, from the fact that different E_ref values give different curves for E_p and I_p, we find that not all transitions are 0–0 vibronic, and thus, again, the Λ values may not be well defined. Nevertheless, for reference, the 0–0 Lambda values range from Λ = 1414.21 for a₁ = 0 to Λ = 0.090 for a₁ = 1.65, E_ref = 12.63 and Λ = 1.25 for a₁ = 1.65, E_ref = 16.84. If these Λ values are relevant, then as in Fig. 2, we are largely operating in the reversible and quasi-reversible regimes.

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 $Ũ_{A}^{μ}$ and $Ũ_{B}^{μ^{'}}$ 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 $\bar{g}$ _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 $\bar{g}$ _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.

FIG. 7.

View large Download slide

1D Linear sweep voltammograms at different values of $\bar{g}$ _p and σ, the mean and variance, respectively, of the lognormal distribution from which g_p values are selected at each time step of the simulation, ω_p = 12.63kT. Simulations are for the harmonic oscillator PESs. The black dotted line shows the case where g_p takes a single value throughout the simulation, equal to the mean of the underlying distribution. The inset of each subplot shows the probability distribution function of g_p values for each parameterization. Due to the stochastic nature of the g_p sampling process, the curves here are averaged over 200 simulations. (a) σ = 0.1, (b) σ = 0.2, and (c) σ = 0.3. As σ approaches zero, the IV curve approaches the single value g case, as expected. However, at larger $\bar{g}$ _p and σ values, we see that the onset of current occurs at lower overpotential. This early onset arises because, as the distribution widens, lower g_p values can be selected with much higher probability, increasing the Franck–Condon factors between the ground vibrational state in 0 and ground vibrational state in 1.

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.³⁹ 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.

FIG. 8.

View large Download slide

Peak current and peak potential as a function of diffusion constant and scan rate for various g_p values for the harmonic oscillator model, ω_p = 12.63kT. [(a) and (c)] Peak current vs diffusion constant and scan rate, respectively. [(b) and (d)] Peak potential vs diffusion constant and scan rate, respectively. The square root dependence of I_p with respect to ν and D applies over a wide range of g_p values, suggesting that this observable is not suitable for determining the presence of proton transfer in an electrochemical ET/PCET experiment. While the behavior of E_p with respect to ν and D does appear to be different depending on g_p, there exists no simple analytical comparison for non-Nernstian ET and thus no means to use Eq. (44) to diagnose whether a given reaction is ET or PCET.

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).

FIG. 9.

View large Download slide

Voltammograms for a quartic potential PCET simulation for hydrogen and deuterium for multiple values of a₁. Voltammograms for hydrogen and deuterium, with L = 0.204, E_ref = 12.63kT, a₄ = 0.005, and a₂ = 0.05. The traces, going for left to right, are for increasing values of a₁. Blue traces correspond to hydrogen and red traces correspond to deuterium. Note that, for a₁ = 0, hydrogen and deuterium yield equivalent voltammetry curves, as expected. While I_p shows almost identical behavior regardless of a₁ and is not drastically different between hydrogen and deuterium, E_p is considerably larger for the hydrogen case compared to deuterium, suggesting that there are significant isotope effects on E_p. The difference in E_p between hydrogen and deuterium is independent of the scan rate and diffusion constant (not shown). Note that the peak position can shift by up to almost 150 mV. The Λ values for these results range from Λ = 1414.21 for a₁ = 0 to Λ = 0.09 for a₁ = 0.006 for hydrogen and Λ = 0.000 65 for a₁ = 0.006 for deuterium. For reference, the Λ values for a₁ = 0.002 are Λ = 106.06 for hydrogen and Λ = 16.24 for deuterium.

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₁ values from the quartic potential model, as well as a plot of $E_{p}^{P C E T (D)} - E_{p}^{P C E T (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 $U_{0}^{μ μ^{'}}$ 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); $ρ (t) = ρ_{b}^{e q} \otimes ρ_{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.

FIG. 10.

View large Download slide

1D Linear sweep voltammograms showing agreement between 1D and 2D simulations for a case with two vibrational states. ω_p = 14.74kT, D = 4.75 · 10⁻⁴, ν = 10⁻¹⁶, ω = 0.000 25, Γ = 10⁻⁵. All parameters are in atomic units. The 2D simulations are represented by dashed lines and the 1D simulations by solid lines. Since the 2D and 1D simulations are in quantitative agreement, we can simulate using the less computationally intensive 1D model to study these PCET systems.

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.42kT), there is a spacing of over 42kT between the ground and highest excited vibrational state.

FIG. 11.

View large Download slide

1D Linear sweep voltammograms showing agreement between simulations with two and five vibrational states. In both subfigures, the dotted lines correspond to N_μ = 2 and solid semi-transparent lines to N_μ = 5. The black trace shows the same simulation, but with no proton coordinate. (a) ω_p held constant at 12.63kT and (b) g_p held constant at 2.83. Since the N_μ = 2 and N_μ = 5 simulations are in quantitative agreement (the curves are on top of each other in the figure), we surmise that the relevant chemistry occurs between no more than the first two vibrational states, and no additional information is gained or lost by increasing the number of vibrational states.

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 ρ₀ only, as the final result for ρ₁ can easily be solved by analogy. We first start with the analog Eq. (22),

\begin{align} ⟨μ| {\hat{\hat{L}}}^{v i b} {\hat{ρ}}_{0}^{v i b} |ν⟩ \\ = - \sum_{k, μ^{'}, ν^{'}, ξ} \frac{1}{ℏ^{2}} \int_{0}^{\infty} d τ (e^{- \frac{i ϵ_{k} τ}{ℏ}} f (ϵ_{k}) ⟨μ| ρ_{0} e^{- \frac{i H_{0} τ}{ℏ}} |ξ⟩ ⟨ξ| W_{k} |μ^{'}⟩) \\ + (⟨μ^{'}| e^{\frac{i H_{1} τ}{ℏ}} W_{k} |ν⟩) \\ - e^{- \frac{i ϵ_{k} τ}{ℏ}} (1 - f (ϵ_{k})) ⟨μ| e^{- \frac{i H_{0} τ}{ℏ}} W_{k} e^{\frac{i H_{1} τ}{ℏ}} |μ^{'}⟩ ⟨μ^{'}| ρ_{1} |ν^{'}⟩ ⟨ν^{'}| W_{k} |ν⟩ \\ - e^{\frac{i ϵ_{k} τ}{ℏ}} (1 - f (ϵ_{k})) ⟨μ| W_{k} |μ^{'}⟩ ⟨μ^{'}| ρ_{1} e^{- \frac{i H_{1} τ}{ℏ}} |ν^{'}⟩ ⟨ν^{'}| W_{k} e^{\frac{i H_{0} τ}{ℏ}} |ν⟩ \\ + (e^{\frac{i ϵ_{k} τ}{ℏ}} f (ϵ_{k}) ⟨μ| W_{k} |μ^{'}⟩ ⟨μ^{'}| e^{- \frac{i H_{1} τ}{ℏ}} W_{k} |ξ⟩ ⟨ξ| e^{\frac{i H_{0} τ}{ℏ}} ρ_{0} |ν⟩) . \end{align}

(C1)

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

\begin{align} \frac{\partial ρ_{0}^{μ ν}}{\partial t} & = \frac{1}{2} ⟨μ| ({Ĥ_{0}^{v i b}, {\hat{ρ}}_{0}^{v i b}} - {{\hat{ρ}}_{0}^{v i b}, Ĥ_{0}^{v i b}}) |ν⟩ - \frac{i}{ℏ} (U_{0}^{μ} - U_{0}^{ν}) ρ_{0}^{μ ν} \\ - \frac{π}{ℏ} \sum_{k, μ^{'}, ξ} f (U_{1}^{μ^{'}} - U_{0}^{ξ}) ⟨μ| W_{k} |μ^{'}⟩ ⟨μ^{'}| W_{k} |ξ⟩ ρ_{0}^{ξ ν} \\ - \frac{π}{ℏ} \sum_{k, μ^{'}, ξ} ρ_{0}^{μ ξ} f (U_{1}^{μ^{'}} - U_{0}^{ξ}) ⟨ξ| W_{k} |μ^{'}⟩ ⟨μ^{'}| W_{k} |ν⟩ \\ + \frac{π}{ℏ} \sum_{k, μ^{'}, ν^{'}} (1 - f (U_{1}^{ν^{'}} - U_{0}^{ν})) ⟨μ| W_{k} |μ^{'}⟩ ⟨ν^{'}| W_{k} |ν⟩ ρ_{1}^{μ^{'} ν^{'}} \\ + \frac{π}{ℏ} \sum_{k, μ^{'}, ν^{'}} ρ_{1}^{μ^{'} ν^{'}} (1 - f (U_{1}^{μ^{'}} - U_{0}^{μ})) ⟨μ| W_{k} |μ^{'}⟩ ⟨ν^{'}| W_{k} |ν⟩ . \end{align}

(C2)

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

Γ_{ν ν^{'}}^{μ μ^{'}} = 2 π \sum_{k} ⟨μ| W_{k} |μ^{'}⟩ ⟨ν| W_{k} |ν^{'}⟩ δ (ϵ - ϵ_{k}),

(C3)

which when is substituted into Eq. (C2) yields

\begin{align} \frac{\partial ρ_{0}^{μ ν}}{\partial t} & = \frac{1}{2} ⟨μ| ({Ĥ_{0}^{v i b}, {\hat{ρ}}_{0}^{v i b}} - {{\hat{ρ}}_{0}^{v i b}, Ĥ_{0}^{v i b}}) |ν⟩ - \frac{i}{ℏ} (U_{0}^{μ} - U_{0}^{ν}) ρ_{0}^{μ ν} \\ - \frac{1}{2 ℏ} \sum_{μ^{'}, ξ} Γ_{μ^{'} ξ}^{μ μ^{'}} f (U_{1}^{μ^{'}} - U_{0}^{ξ}) ρ_{0}^{ξ ν} - \frac{1}{2 ℏ} \sum_{μ^{'}, ξ} Γ_{μ^{'} ν}^{ξ μ^{'}} ρ_{0}^{μ ξ} f (U_{1}^{μ^{'}} - U_{0}^{ξ}) \\ + \frac{1}{2 ℏ} \sum_{μ^{'}, ν^{'}} Γ_{ν^{'} ν}^{μ μ^{'}} (1 - f (U_{1}^{ν^{'}} - U_{0}^{ν})) ρ_{1}^{μ^{'} ν^{'}} \\ + \frac{1}{2 ℏ} \sum_{μ^{'}, ν^{'}} Γ_{ν^{'} ν}^{μ μ^{'}} ρ_{1}^{μ^{'} ν^{'}} (1 - f (U_{1}^{μ^{'}} - U_{0}^{μ})) . \end{align}

(C4)

Finally, taking the secular approximation, we arrive at

\begin{align} \frac{\partial ρ_{0}^{μ μ}}{\partial t} & = ⟨μ| ({Ĥ_{0}^{v i b}, {\hat{ρ}}_{0}^{v i b}}) |μ⟩ - \frac{1}{ℏ} \sum_{μ^{'}} Γ_{μ^{'} μ}^{μ μ^{'}} f (U_{1}^{μ^{'}} - U_{0}^{μ}) ρ_{0}^{μ μ} \\ + \frac{1}{ℏ} \sum_{μ^{'}} Γ_{μ^{'} μ}^{μ μ^{'}} (1 - f (U_{1}^{μ^{'}} - U_{0}^{μ})) ρ_{1}^{μ^{'} μ^{'}}, \end{align}

(C5)

which is similar to Eq. (26), except that $Γ {|F_{μ \to μ^{'}}|}^{2}$ has been replaced with the more general $Γ_{μ^{'} μ}^{μ μ^{'}}$ ⁠.

REFERENCES

1.

G.

Herzberg

,

Molecular Spectra and Molecular Structure

(

Van Nostrand

,

New York

,

1950

).

Google Scholar

2.

H.

Sumi

and

R. A.

Marcus

,

J. Chem. Phys.

84

,

4894

(

1986

).

https://doi.org/10.1063/1.449978

Crossref

3.

J.

Lappe

and

R. J.

Cave

,

J. Phys. Chem. A

104

,

2294

(

2000

).

https://doi.org/10.1021/jp992518z

Crossref

4.

X.

Sun

and

E.

Geva

,

J. Phys. Chem. A

120

,

2976

(

2016

).

https://doi.org/10.1021/acs.jpca.5b08280

Crossref

PubMed

5.

E.

Buhks

,

M.

Bixon

,

J.

Jortner

, and

G.

Navon

,

J. Phys. Chem.

85

,

3759

(

1981

).

https://doi.org/10.1021/j150625a010

Crossref

6.

J. A. V.

Butler

,

Proc. R. Soc. London, Ser. A

157

,

423

(

1936

).

https://doi.org/10.1098/rspa.1936.0205

Crossref

7.

D. M.

Newns

,

Phys. Rev.

178

,

1123

(

1969

).

https://doi.org/10.1103/physrev.178.1123

Crossref

8.

W.

Schmickler

,

J. Electroanal. Chem.

204

,

31

(

1986

).

https://doi.org/10.1016/0022-0728(86)80505-8

Crossref

9.

A.

Soudackov

and

S.

Hammes-Schiffer

,

J. Chem. Phys.

111

,

4672

(

1999

).

https://doi.org/10.1063/1.479229

Crossref

10.

A.

Soudackov

and

S.

Hammes-Schiffer

,

J. Chem. Phys.

113

,

2385

(

2000

).

https://doi.org/10.1063/1.482053

Crossref

11.

A.

Soudackov

,

E.

Hatcher

, and

S.

Hammes-Schiffer

,

J. Chem. Phys.

122

,

014505

(

2005

).

https://doi.org/10.1063/1.1814635

Crossref

12.

J.

Grimminger

,

S.

Bartenschlager

, and

W.

Schmickler

,

Chem. Phys. Lett.

416

,

316

(

2005

).

https://doi.org/10.1016/j.cplett.2005.09.084

Crossref

13.

D. R.

Weinberg

,

C. J.

Gagliardi

,

J. F.

Hull

,

C. F.

Murphy

,

C. A.

Kent

,

B. C.

Westlake

,

A.

Paul

,

D. H.

Ess

,

D. G.

McCafferty

, and

T. J.

Meyer

,

Chem. Rev.

112

(

7

),

4016

–

4093

(

2012

).

https://doi.org/10.1021/cr200177j

Google Scholar

Crossref

14.

M. T.

Koper

,

Chem. Sci.

4

,

2710

(

2013

).

https://doi.org/10.1039/c3sc50205h

Crossref

15.

C.

Costentin

and

J.-M.

Savéant

,

J. Am. Chem. Soc.

139

,

8245

(

2017

).

https://doi.org/10.1021/jacs.7b02879

Crossref

PubMed

16.

P.

Goyal

and

S.

Hammes-Schiffer

,

J. Phys. Chem. Lett.

6

,

3515

(

2015

).

https://doi.org/10.1021/acs.jpclett.5b01475

Crossref

PubMed

17.

P.

Goyal

,

C. A.

Schwerdtfeger

,

A. V.

Soudackov

, and

S.

Hammes-Schiffer

,

J. Phys. Chem. B

120

,

2407

(

2016

).

https://doi.org/10.1021/acs.jpcb.5b12015

Crossref

PubMed

18.

K.

Song

and

Q.

Shi

,

J. Chem. Phys.

146

,

184108

(

2017

).

https://doi.org/10.1063/1.4982928

Crossref

19.

P.

Goyal

and

S.

Hammes-Schiffer

,

Proc. Natl. Acad. Sci. U. S. A.

114

,

1480

(

2017

).

https://doi.org/10.1073/pnas.1621393114

Crossref

PubMed

20.

A.

Mandal

,

F. A.

Shakib

, and

P.

Huo

,

J. Chem. Phys.

148

,

244102

(

2018

).

https://doi.org/10.1063/1.5030634

Crossref

PubMed

21.

S.

Ghosh

,

J.

Castillo-Lora

,

A. V.

Soudackov

,

J. M.

Mayer

, and

S.

Hammes-Schiffer

,

Nano Lett.

17

,

5762

(

2017

).

https://doi.org/10.1021/acs.nanolett.7b02642

Crossref

PubMed

22.

C.

Costentin

,

M.

Robert

, and

J.-M.

Savéant

,

J. Am. Chem. Soc.

128

(

14

),

4552

–

4553

(

2006

).

https://doi.org/10.1021/ja060527x

Crossref

PubMed

23.

M. T.

Huynh

,

S. J.

Mora

,

M.

Villalba

,

M. E.

Tejeda-Ferrari

,

P. A.

Liddell

,

B. R.

Cherry

,

A.-L.

Teillout

,

C. W.

MacHan

,

C. P.

Kubiak

,

D.

Gust

,

T. A.

Moore

,

S.

Hammes-Schiffer

, and

A. L.

Moore

,

ACS Cent. Sci.

3

,

372

(

2017

).

https://doi.org/10.1021/acscentsci.7b00125

Crossref

PubMed

24.

W.

Dou

and

J. E.

Subotnik

,

J. Chem. Phys.

145

,

054102

(

2016

).

https://doi.org/10.1063/1.4959604

Crossref

PubMed

25.

A. J.

Coffman

,

A. K.

Harshan

,

S.

Hammes-Schiffer

, and

J. E.

Subotnik

,

J. Phys. Chem. C

123

,

13304

(

2019

).

https://doi.org/10.1021/acs.jpcc.9b02068

Crossref

26.

H_sys can also be expressed in second quantized notation as $H_{s y s} = \sum_{μ} U_{0}^{μ} (1 - d^{†} d) a_{μ}^{†} a_{μ} + \sum_{μ^{'}} U_{1}^{μ^{'}} d^{†} d a_{μ^{'}}^{†} a_{μ^{'}} + H_{p h}$ ⁠, where $a_{μ} (a_{μ}^{†})$ are annihilation (creation) operators for the proton vibrational state μ and $1$ is the identity operator.

27.

K.

Blum

,

Density Matrix Theory and Applications

, 3rd ed. (

Springer-Verlag Berlin Heidelberg

,

2012

).

Google Scholar

Crossref

28.

For a previous application of the QME to photo-induced PCET, see Refs. 47 and 48.

29.

W.

Dou

and

J. E.

Subotnik

,

J. Chem. Phys.

144

,

024116

(

2016

).

https://doi.org/10.1063/1.4939734

Crossref

PubMed

30.

R.

Kapral

,

J. Chem. Phys.

110

,

8919

(

1999

).

https://doi.org/10.1063/1.478811

Crossref

31.

Y. J.

Jung

and

J.

Cao

,

J. Chem. Phys.

117

,

3822

(

2002

).

https://doi.org/10.1063/1.1491241

Crossref

32.

I.

Navrotskaya

,

A. V.

Soudackov

, and

S.

Hammes-Schiffer

,

J. Chem. Phys.

128

,

244712

(

2008

).

https://doi.org/10.1063/1.2940203

Crossref

PubMed

33.

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.

34.

C.

Venkataraman

,

A. V.

Soudackov

, and

S.

Hammes-Schiffer

,

J. Phys. Chem. C

112

,

12386

(

2008

).

https://doi.org/10.1021/jp802171y

Crossref

35.

Z.

Zhou

,

H.-T.

Chen

,

A.

Nitzan

, and

J. E.

Subotnik

, arXiv:1909.11695 (

2019

).

36.

We note that recent work by Matyushov and Newton⁴⁹ 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 λ.

37.

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

38.

H.

Matsuda

and

Y.

Ayabe

,

Z. Electrochem

59

,

494

(

1955

).

39.

A. J.

Bard

and

L. R.

Faulkner

,

Electrochemical Methods: Fundamentals and Applications

(

Wiley

,

New York

,

2001

).

Google Scholar

40.

J.

Koch

,

F.

Von Oppen

,

Y.

Oreg

, and

E.

Sela

,

Phys. Rev. B

70

,

195107

(

2004

).

https://doi.org/10.1103/physrevb.70.195107

Crossref

41.

J.

Koch

,

M.

Semmelhack

,

F.

Von Oppen

, and

A.

Nitzan

,

Phys. Rev. B

73

,

155306

(

2006

).

https://doi.org/10.1103/physrevb.73.155306

Crossref

42.

J. P.

Dahl

and

M.

Springborg

,

J. Chem. Phys.

88

,

4535

(

1988

).

https://doi.org/10.1063/1.453761

Crossref

43.

A.

Migliore

,

N. F.

Polizzi

,

M. J.

Therien

, and

D. N.

Beratan

,

Chem. Rev.

114

,

3381

(

2014

).

https://doi.org/10.1021/cr4006654

Crossref

PubMed

44.

A. V.

Soudackov

and

S.

Hammes-Schiffer

,

J. Chem. Phys.

143

,

194101

(

2015

).

https://doi.org/10.1063/1.4935045

Crossref

PubMed

45.

Z. K.

Goldsmith

,

A. V.

Soudackov

, and

S.

Hammes-Schiffer

,

Faraday Discuss.

216

,

363

(

2019

).

https://doi.org/10.1039/c8fd00240a

Crossref

PubMed

46.

a₄ and a₂ can also change the observed IV behavior but do not have as simple of physical intuition as the a₁ 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₁ dependence, as this illustrates qualitatively how the quartic oscillator shows similar current vs potential trends as the other oscillators examined.

47.

C.

Venkataraman

,

A. V.

Soudackov

, and

S.

Hammes-Schiffer

,

J. Chem. Phys.

131

,

154502

(

2009

).

https://doi.org/10.1063/1.3249964

Crossref

PubMed

48.

C.

Venkataraman

,

A. V.

Soudackov

, and

S.

Hammes-Schiffer

,

J. Phys. Chem. C

114

,

487

(

2010

).

https://doi.org/10.1021/jp908798n

Crossref

49.

D. V.

Matyushov

and

M. D.

Newton

,

J. Chem. Phys.

147

,

194506

(

2017

).

https://doi.org/10.1063/1.5003022

Crossref

PubMed

2020

Author(s)

Modeling voltammetry curves for proton coupled electron transfer: The importance of nuclear quantum effects

I. INTRODUCTION

II. THEORY AND METHODS

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

1. Anderson Holstein (AH) Hamiltonian for PCET

2. Redfield perturbation theory in the interaction picture

3. Redfield perturbation theory in the Schrödinger picture

4. A partial Wigner transform

5. A secular approximation

6. Final equations

B. Simulation methods for generating IV curves on a grid

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

2. Simulating voltammetry and measuring current

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

C. Parameters

D. Understanding and interpreting IV curves

III. RESULTS

A. Harmonic oscillator eigenstates

B. Morse potential

C. Quartic potential

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

IV. DISCUSSION

A. Dependence on diffusion constant and scan rate

B. Hydrogen/deuterium isotope effects in a quartic potential

V. CONCLUSION

ACKNOWLEDGMENTS

APPENDIX A: EQUIVALENCE OF 2D AND 1D SIMULATIONS

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

APPENDIX C: W_k DEPENDENCE ON PROTON COORDINATE

REFERENCES

Citing articles via

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

Modeling voltammetry curves for proton coupled electron transfer: The importance of nuclear quantum effects

I. INTRODUCTION

II. THEORY AND METHODS

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

1. Anderson Holstein (AH) Hamiltonian for PCET

2. Redfield perturbation theory in the interaction picture

3. Redfield perturbation theory in the Schrödinger picture

4. A partial Wigner transform

5. A secular approximation

6. Final equations

B. Simulation methods for generating IV curves on a grid

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

2. Simulating voltammetry and measuring current

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

C. Parameters

D. Understanding and interpreting IV curves

III. RESULTS

A. Harmonic oscillator eigenstates

B. Morse potential

C. Quartic potential

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

IV. DISCUSSION

A. Dependence on diffusion constant and scan rate

B. Hydrogen/deuterium isotope effects in a quartic potential

V. CONCLUSION

ACKNOWLEDGMENTS

APPENDIX A: EQUIVALENCE OF 2D AND 1D SIMULATIONS

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

APPENDIX C: Wk DEPENDENCE ON PROTON COORDINATE

REFERENCES

Citing articles via

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

This Feature Is Available To Subscribers Only

APPENDIX C: W_k DEPENDENCE ON PROTON COORDINATE