A Grid-free Approach for Simulating Sweep and Cyclic Voltammetry
AA Grid-free Approach for Simulating Sweep and Cyclic Voltammetry.
Alec J. Coffman, Jianfeng Lu, and Joseph E. Subotnik Department of Chemistry, University of Pennsylvania, Philadelphia,Pennsylvania 19104, USA Departments of Mathematics, Physics, and Chemistry, Duke University, Durham,North Carolina 27708, USA (Dated: 10 February 2021)
We present a new computational approach to simulate linear sweep and cyclic voltam-metry experiments that does not require a discretized grid in space to quantify dif-fusion. By using a Green’s function solution coupled to a standard implicit ordinarydifferential equation solver, we are able to simulate current and redox species concen-trations using only a small grid in time . As a result, where benchmarking is possible,we find that the current method is faster (and quantitatively identical) to establishedtechniques. The present algorithm should help open the door to studying adsorptioneffects in inner sphere electrochemistry.1 a r X i v : . [ phy s i c s . c h e m - ph ] F e b . INTRODUCTIONA. Models for Outer-Sphere Electrochemical electron transfer (ET) andVoltammetry A common experimental technique throughout the field of electrochemistry is sweepvoltammetry, where the electrostatic potential of a metal electrode (relative to an electro-chemical cell containing an electroactive redox species) is changed with time. The resultantcurrent is usually plotted as a function of this change in driving force . When the potentialis ramped in a linear fashion, E ( t ) = E + νt , (where ν is the scan rate), the resultingplot is known as a linear sweep voltammetric (LSV) curve. When the scan is carried outin both directions (oxidation and reduction), it is known as cyclic voltammetry (CV). Thecurrent in such experiments arise from both (i) the change in the electron transfer (ET)kinetics at the electrode (due to changing the driving force for the ET reaction), and also(ii) diffusional effects which create a concentration gradient of the redox active species insolution. This novel combination of effects driving the current vs. potential profile makeslinear sweep voltammetry an ideal experiment for inferring unknown chemical and physicalparameters of interest in the electrochemical system. In particular, nowadays, it is routineto use LSV to infer the diffusion constant of the redox active species or the standard rateconstant for an ET reaction, k (the rate constant where the forward and backward rate forthe redox reactions are equal) .In order to properly analyze an LSV experiment, one requires a mathematical modelwith a reasonable description of both the diffusional and ET effects. Consider the mostbasic electrochemical reaction, the reduction of a redox species A to its reduced form B ,A + e – k f −− (cid:42)(cid:41) −− k b B .where species A exists in bulk solution but, before a reaction, must be brought to the metalsurface. For most simulations, diffusional effects are modeled through a standard parabolic,second order partial differential equation (PDE) , ∂c A ∂t = D A ∂ c A ∂x ∂c B ∂t = D B ∂ c B ∂x . (1)Thereafter (and see Eq. 8 below), ET effects are incorporated through either a source term inthe PDE or electrode boundary conditions. This approach has wide spread applicability and2an be adadpted to take into account electrical double layer (EDL) effects (which are alsocalled Frumkin corrections to the rate constant that account for uncompensated resistancedue to field drop at the point/plane of ET) , surface adsorption of a species pre-ET ,coupled homogenous chemical reactions , and many other effects as well. Some of therelevant PDEs can be solved analytically if the ET kinetics are considered to be infinitely fast(i.e. when ET is taken to be instantaneous and dependent upon only the thermodynamics ofthe redox couple at a given external potential); this is known as the Nernstian limit (again,see below). Realistically, however, most ET kinetics are not infinitely fast and the Nernstianlimit is an idealization. In practice, in order to interpret his/her data, the experimentalistwill usually need to model voltammetry experiments explicitly, i.e. one requires a computerto solve the necessary PDEs (e.g., in Eq. 8) that incorporate (at minimum) both diffusionand a time-dependent ET rate constant.Finally, for completeness, a few words are appropriate as far as how most electrochemicalsimulations model electron transfer if one wishes to go beyond the Nernstian limit. Ingeneral, there are two common choices for the rate constant: Butler-Volmer (BV) andMarcus-Hush-Chidsey (MHC).1. The most common choice for the ET rate constant is the Butler-Volmer (BV) ex-pression. BV dynamics are characterized by two parameters: the charge transfercoefficient, α , and the standard prefactor k . The charge transfer coefficient’s physicalinterpretation has historically been a matter of some debate , but mathematically isdefined as the change in the log of the current with respect to change in the drivingforce (in unitless form): α = (cid:12)(cid:12) RTF d ln IdE (cid:12)(cid:12) . The total BV forward ( k f ) and backward ( k b )rate constants take the form k f = k e − α FRT ( E − E (cid:48) ) k b = k e (1 − α ) FRT ( E − E (cid:48) ) , (2)where E (cid:48) is the formal potential (which is dependent on the experimental setup), F is the Faraday constant, and R is the universal gas constant.2. The less common form for the ET rate constant is the MHC expression. Like BV, MHCcan also be characterized by two parameters, the prefactor k , and the reorganizationenergy λ . While the prefactor k has the same interpretation as in BV theory, the3 G b FIG. 1. The standard harmonic potential energy surfaces over a general solvent reorganizationcoordinate within the Marcus picture. The red curve represents the charge located on a redoxspecies molecule, while the blue curves represent the charge being located on one of the continuumof possible electronic states that are accessible in the metal electrode. The free energy differencebetween the charge being located on the molecule versus the metallic state at the fermi level ofthe metal electrode is given by ∆ G , while the reorganization energy ( λ ) is the energy required tochange the relaxed conformation of the redox active species from a local solvent environment ofreactant state to the local solvent environment in the product state (all while keeping the chargefixed in the product state). The barrier for transferring an electron from the product to reactantstate is also shown, labeled by ∆ G ‡ b . reorganization energy λ is defined as the free energy change associated with relaxing thelocal solvent environment from one configuration (equilibrated to the reactant state)to another configuration (equilibrated to the product state), all the while keeping themolecule’s electronic state fixed as the product.4 For homogenous ET in solution, at high temperature, MHC theory is equivalentto transition state theory, such that: k = k e − ∆ G ‡ kBT (3)where ∆ G ‡ is the barrier height associated with the transition state, while k B and T have their usual definitions as the Boltzmann constant and temperature,respectively . While the exact form of k can be complicated and depend onthe nature of the reaction (inner or outer sphere) , such a distinction will notbe relevant for our present purposes. In practice, if one makes a parabolic freeenergy assumption, the Marcus rate takes on the famous form: k = k e ( λ +∆ G )24 kBT , (4)where ∆ G is the free energy difference between the electron being located on theacceptor versus the donor. • For a heterogeneous electrochemical ET, the Marcus rate constant must be mod-ified to account for the fact that the product is no longer a chemical species insolution, but rather a metal electrode with a continuum of electronic states, whichobeys an occupancy dictated by a Fermi distribution at a given external poten-tial. By integrating the rate constants over the distribution of electronic statesin the electrode, one obtains the so-called MHC rate constant for electrochemicalET, k f = k (cid:90) ∞−∞ f ( (cid:15) ) e − ∆ G ‡ f ( (cid:15) ) RT d(cid:15)k b = k (cid:90) ∞−∞ f ( − (cid:15) ) e − ∆ G ‡ b ( (cid:15) ) RT d(cid:15), (5)Here f ( (cid:15) ) = e β(cid:15) is the fermi function and β = k B T . B. Adsorptive Effects in Linear Sweep Voltammetry-A First Step TowardsInner-Sphere ET
The above description of ET in LSV experiments is usually sufficient for ET that occursin solution between redox species. However, if ET occurs between redox species that are5dsorbed to the surface of the electrode, many new factors are introduced that complicatethe mathematics. These factors include (i) possible changes in the driving force for ET dueto the adsorbed species being located within the electrical double layer (EDL); (ii) possiblydifferent reaction mechanisms that can affect both the kinetics and thermodynamics, i.e.inner vs. outer sphere, adiabatic vs. nonadiabatic mechanisms; (iii) the introduction of anadditional timescale for mass transfer, i.e. one must now account not only for a species inthe bulk approaching the surface, but also for the rate of adsorption and desorption of theoxidized and reduced species (as well as a diffusion times moving along the surface). Thislist of important focuses at a surface is by no means exhaustive, but these three effects arelikely the most important.Historically, the numerical simulation of adsorption effects in voltammetry experimentshas been challenging because adsorption introduces nonlinearity into the differential equationin Eq. 1; see Eq. 8 below. In order to conceptually separate nonlinear adsorption fromthe actual ET event, electrochemists have developed a few different model isotherms thatquantify what concentration of a species A is expected to be adsorbed on a given surface(Γ A ) given a bulk concentration in solution ( c A ). The most common adsorption isothermis the Langmuir isotherm , which assumes that, at high enough solution concentrations,there is a saturating coverage of the electrode (denoted Γ s ). The Langmuir isotherm furtherassumes that there are no interactions between the species on the electrode surface andthat there is no surface heterogeneity. Mathematically, for a single species A , the Langmuirisotherm is given by: Γ A Γ s − Γ A = β A c A ⇐⇒ Γ A Γ s = β A c A β A c A , (6)where β A is the equilibrium constant for the adsorption/desorption of A , β A = k Aads k Ades . Forthe case of multiple species (for example, A and B , where B is the reduced form of A ), thesurface concentration of species i ( i = A, B ) can be expressed asΓ i = Γ s β i c i β i c i + (cid:80) j (cid:54) = i β j c j (7)More complicated isotherms are sometimes used in the literature, such as the Frumkinisotherm (which allows for interactions between adsorbed species and the surface which canbe coverage/concentration dependent). Below, we will present data from only the Langmuirisotherm (though the method itself will be quite general).6 I. ALGORITHMA. Discretized Method
Before presenting our grid-free approach, we quickly review the standard approach forelectrochemical voltammetry in the absence of adsorption. Assuming that all electron trans-fer occurs only at the boundary, the differential equation for the time evolution of the con-centration of each species is: ∂c A ∂t = D A ∂ c A ∂x − ( k f c A − k b c B ) δ ( x ) ∂c B ∂t = D B ∂ c B ∂x + ( k f c A − k b c B ) δ ( x ) , (8)Here, c A and c B are the normalized concentrations of species A and B , respectively, and D A and D B are the diffusion constants for species A and B . The delta function source term isused to enforce that ET occurs only at the solution-electrode interface, defined as x = 0.The equations above are typically solved on a discretized one-dimensional grid. In otherwords, for a discrete set of points { x , x , ..., x N } with spacing ∆ x and a time step ∆ t , onesolves for c ( j ∆ x, k ∆ t ) in a big loop. In principle, it is usually trivial to solve the differentialequation in Eq. 8 using the Euler method with a small enough time step. Unfortunately,however, the ET source terms at the metal electrode (i.e. those terms proportional to δ ( x ))usually lead to a myriad of numerical problems for any realistic time discretization (thatwould certainly dissuade most undergraduates). Moreover, we cannot emphasize enoughthat electrochemistry occurs over a vast array of time scales. While electron transfer mayoccur on the time scale of picoseconds, an electrochemical sweep over 1 V with a scan rateof 10mV/s occurs over 100 seconds; one must propagate Eq. 8 over an immensely long timerange.The most common solution to such numerical difficulties is to transform the source termsin the above equation into boundary conditions, utilizing a no-flux boundary condition atthe electrode surface; for details see Ref. 18. Thereafter, one must still invoke an implicitODE solver, as the equations are stiff . Alternatively, because the operator on the righthand side of Eq. 8 is linear, the simplest approach is to just diagonalize the corresponding(non-Hermitian) operator for diffusion and ET. While diagonalization of a non-Hermitianoperator is usually unstable, it is well known that the non-Hermitian diffusion operator(in one dimension) can be transformed into a Hermitian operator . As a result, if one is7repared to work with a grid of points in one-spatial dimension, and one is not concernedwith adsoprtion, there are a few numerically stable techniques available as far calculatingLSV curves. Our recent work has shown that one can use such a brute force approachto capture more complicated electrochemical scenarios, including proton coupled electontransfer . B. Grid-Free Method
Unfortunately, the methods above face great difficulties when nonlinear adsorption anddesorption are considered. On the one hand, the presence of non-linear Langmuir isothermslead to an even stronger stiffening of the coupled set of PDEs, such that using an implicitdifferntial equation solver to solve for the dynamics on a grid often fails. On the other hand,diagonalization is obviously impossible with nonlinear equations of motion. As a result,modeling adsorption is daunting for most parameter sets. The most successful algorithm todate comes from Compton et al who have used a spatial grid in a recent article to studyphase transitions of adsorbed electroactive species .The goal of the present letter is to present a different solution to this problem, one thatcompletely removes the necessity for a grid in space. In so doing, we will vastly reduce thenumber of variables that must be propagated in time and, rather than spend time solving theheat equation in solution (and calculating the free solution concentration c A ), we will focusour attention exclusively at the surface where the electron transfer occurs (and calculate thesurface concentration Γ A ).
1. The relevant PDEs
Before outlining our new approach, let us define the differential equations for couplingdiffusion, ET, and adsorption that we wish to solve. The PDEs for the solution concentra-tions ( c A , c B ) take the standard form (following Eq. 8 above), with the addition of adsorp-tion/desorption rates, ν ads / ν des at the electrode surface ( x = 0): ∂c A ∂t = D A ∂ c A ∂x − [( k solf ( t ) c A − k solb ( t ) c B ) + ν ads,A − ν des,A ] δ ( x ) (9a) ∂c B ∂t = D B ∂ c B ∂x + [( k solf ( t ) c A − k solb ( t ) c B ) − ν ads,B + ν des,B ] δ ( x ) . (9b)8he ODEs for the surface concentrations, Γ A and Γ B involve both electron transfer as wellas surface adsorption/desoprtion: ∂ Γ A ∂t = ν ads,A ( t ) − ν des,A ( t ) − ( k adsf ( t )Γ A ( t ) − k adsb ( t )Γ B ( t )) (10a) ∂ Γ B ∂t = ν ads,B ( t ) − ν des,B ( t ) + ( k adsf ( t )Γ A ( t ) − k adsb ( t )Γ B ( t )) . (10b)The choice of isotherm will dictate the exact form of ν iads /ν ides ; e.g. for the case of a Langmuirisotherm, we choose ν ads,A ( t )= k Aads c A ( t, x = 0)[Γ s − (Γ A ( t ) + Γ B ( t ))] (11a) ν des,A ( t )= k Ades Γ A ( t ) (11b) ν ads,B ( t )= k Bads c B ( t, x = 0)[Γ s − (Γ A ( t ) + Γ B ( t ))] (11c) ν des,B ( t )= k Bdes Γ B ( t ) . (11d) k adsf /k solf ( k adsb /k solb ) can follow different ET rate laws and are parameterized in different ways(different α/λ, k , E , etc.) to reflect the correct underlying physics.Lastly, when calculating the current, note that in principle there will be separate contri-butions from both solution and adsorbed ET, I ( t ) = nF A [( k solf ( t ) c A ( t, x = 0) − k solb ( t ) c B ( t, x = 0)) + ( k adsf ( t )Γ A − k adsb ( t )Γ B )] , (12)where n is the number of electrons, A the surface area of the electrode, and F is theFaraday constant. Formally, the current depends only on four variables at a single time: c A ( x = 0 , t ) , c B ( x = 0 , t ) , Γ A ( t ) , and Γ B ( t ).
2. A Green’s Function Approach
Given that ν ads,i /ν des,i and I ( t ) depend only on c A ( x = 0) , c B ( x = 0) , Γ A , and Γ B (andspecifically not on any of the remaining c A ( x (cid:54) = 0) , c B ( x (cid:54) = 0)), our goal is to find a maximallyefficient algorithm that tracks only these four values (without tracking the remaining solutionconcentrations away from the electrode surface). Here, for simplicity, we will assume thatspecies A and B have equal diffusivities ( D A = D B ); all of the following results can besimply modified to accommodate unequal diffusion constants, if needed.To avoid using a grid in space, we invoke Duhamel’s principle to recast our coupleddifferential equations modeling c A and c B as integral equations. For an inhomogeneous9iffusion equation of the form Eq. 9, with initial and boundary conditions c B ( t = 0 , x ) = c B ( t, x = ±∞ ) = 0, c A ( t = 0 , x ) = c A ( t, x = ±∞ ) = 1, one can formally invert theseequations of motion to solve for the solution concentration in integral form (we will show c B , but the derivation is identical for c A ): c B ( t, x ) = (cid:90) t (cid:90) x G ( x − y, t − s )[( k solf ( s ) c A ( s, y ) − k solb ( s ) c B ( s, y )) − ν ads,B + ν des,B ] δ ( y ) ds dy = (cid:90) t (cid:90) x G ( x − y, t − s )[ b ( s, y ) + a ( s ) c B ( s, y )] δ ( y ) ds dy = (cid:90) t G ( x, t − s )[ b ( s, y ) + a ( s ) c B ( s, ds (13)where G ( x − y, t − s ) = 1 (cid:112) D B π ( t − s ) e − | x − y | t − s ) (14)is one half of the heat equation Greens function in one dimension . Note that a factor of 2appears here because the original diffusion equation spans only the half-space x ∈ [0 , ∞ ).In Eq. 13, the expression for b ( s, y ) and a ( s ) are obtained directly from reorganization ofEq. 9; for a Langmuir isotherm, these expressions are b ( s, y ) = k solf ( s ) c A ( s, y ) + k Bdes Γ B ( s ) a ( s ) = − ( k solb ( s ) + k Bads [Γ s − (Γ A ( s ) + Γ B ( s ))]) . (15)Since we are only concerned with the solution concentration at x = 0, the equation we wishto numerically solve is given by c B ( t,
0) = (cid:90) t (cid:112) D B π ( t − s ) [ b ( s,
0) + a ( s ) c B ( s, ds. (16)We face two challenges to employing Eq. 16 to measure electrochemical currents:1. b ( s,
0) and a ( s ) depend on Γ A , Γ B , and c A explicitly, while Γ A and Γ B depend explicitlyon c A and c B .2. There is a singularity in Eq. 16 at s = t , which must be carefully accounted for dueto the importance of accurate integral evaluation at this point.10et us begin with the second challenge. To treat the singularity, we cast the integrand inEq. 16 into a general form of f ( s ) √ t − s , where f ( s ) = b ( s,
0) + a ( s ) c B ( s, f ( s ) − f ( t ) √ t − s and f ( t ) √ t − s . The first piece goes to zero at s = t , and the secondpiece can be integrated analytically, removing the singularity. By using nonweighted numer-ical quadrature (according to the trapezoid rule), one can obtain a numerical prescriptionfor calculating c B ( t ) as a function of b ( s, a ( s ), and the prior history of c B (expressed ona grid [ s , s , ..., s N ], with spacing ∆ s ): c B ( t,
0) = ∆ s ( N − (cid:80) i =1 b ( s i , a ( s i ) c B ( s i , − b ( t, √ t − s i + b (0 , − b ( t, √ t ) + 2 b ( t, √ t √ D B π − a ( t )[2 √ t − ∆ s ( √ t + N − (cid:80) i =1 1 t − s i )] . (17)For a derivation, see Appendix I. Note that, in the absence of adsorption, Eq. 17 is entirelysufficient for calculating the current as a function of time/voltage for a sweep voltammetrysimulation – without any need for solving an ODE on a discretized grid in space. Forinstance, suppose that A and B have equal diffusion constants for A and B . In sucha case, c A ( t,
0) + c B ( t,
0) = c A (0 , x = ∞ ) + c B (0 , x = ∞ ) = c bulk , so that c A ( t, x = 0)can be determined by c A ( t,
0) = c bulk − c B ( t, k solf and k solb are analyticexpressions (already needed for evaluation of Eq. 17), one can easily calculate the current attime t using Eq. 12 and 17, and thus challenge .Next, let us consider the first challenge, i.e. how to self consistently calculate concentra-tions when there is adsorption of redox active species. In order to simultaneously solve Eq.10 (a standard ODE) and Eq. 17 (an integro-differential equation) in coupled fashion, wewill apply the procedure below:1. At t = dt , calculate c B (0 , dt ) using Eq. 17 and the values of Γ A , Γ B and c A at theirprior time step values.2. Calculate c A (0 , dt ) using the analogous Eq. 17 for c A , with the value of c B (0 , dt ) cal-culated in step 1 and the values of Γ A and Γ B at their prior time step values.3. Return to step 1-2 until the tolerated convergence is achieved for c A (0 , dt ) and c B (0 , dt ).11. Calculate Γ A ( dt ) and Γ B ( dt ) using Eq. 10 and an implicit ODE solver (in this case,backward differentiation formula of order 5).5. Return to step 1 of this procedure until simultaneous convergence to a predeterminedtolerance is achieved for all four values c A (0 , dt ), c B (0 , dt ), Γ A ( dt ), and Γ B ( dt ) withinone single convergence step.6. Move on to the next time step, repeating steps 1-5 at each time step until the end ofthe simulation is reached.Note that, in step 1) above, for the first guess of Γ A and Γ B at a given time step, extrapolationmethods can be used to obtain a better initial guess than simply using the prior time step;this usually yields much faster convergence. In our experience, linear extrapolation usingtwo prior times (Γ A ( t ) = 2Γ A ( t − δt ) − Γ A ( t − δt )) yields a substantial improvement inconvergence time, while higher order extrapolations do little for accelerating convergence.In Fig. 2 we show a set of illustrative examples of CVs generated with this methodology,both in the case (a) without and (b) with adsorptive effects. From the voltamogramms inFig. 2(a) without adsorption, at low scan rate we see the expected peak separation of 57 mV for a reversible ET reaction in solution, and the peak separation increases slightly as thescan rate increases and approaches the irreversible ET regime. Conversely, in Fig. 2(b),we see no peak separation in the case where ET occurs between surface bound species, asexpected . These results demonstrate that the above algorithm can capture the behavior ofboth coupled mass and electron transfer in solution and at surface bound electrochemicalinterfaces under a sweeping potential.Additionally, to demonstrate the computational benefit of the present algorithm, belowwe show two tables illustrating the computational improvement of this method. Table Ishows computation times for the results shown in Fig. 2, where we see that this presentalgorithm is relatively invariant to scan rate – even with adsorption. This behavior standsin contrast to the typical grid based ODE methods, which can require longer times at smallerscan rates, as more time is needed to propagate dynamics using an implicit ODE solver inorder to cover a given potential range. Lastly, Table II below shows the relative improvementin computation time for the present algorithm as compared to the more traditional approachwhere one solves the spatial grid-based coupled ODEs through discretization of the relevantpropagator over a set of grid points and diagonalization of the diffusion operator . Similar12o the above argument about scan rate for traditional grid based methods, more spatialpoints implies larger computational cost, whereas for the method outlined in present paper,no such discretization is required as there is no added computational cost. All simulationswere performed in Python, using the NumPy and SciPy libraries, on a 2.8 GHz Dual-CoreIntel Core i5 processor with 16 GB of RAM. TABLE I. Computational wall time by method (Python)Scan rate (V/sec) No adsorption Adsorption0.05 3.70 s 22.06 s0.1 4.89s 21.13 s1 3.40 s 20.54 s5 3.67 s 19.87 sTABLE II. Computation time by method (Python)Algorithm wall timeGrid-free 2.34 sDiscretized, N x = 50 9.94 sDiscretized, N x = 100 30.91 sDiscretized, N x = 200 83.02 s III. CONCLUSION
In this article, as an alternative to traditional grid based methods, we have shown anovel integral based approach for simulating electrochemical ET to model ET along withadsorption of redox species. Our approach is in quantitative agreement with grid basedPDE solutions for a single ET experiment; more generally, however, the present algorithmshould allow us to simulate more complicated sweep experiments that involve adsorptionkinetics and beyond. Code implementing this algorithm can be found at https://github.com/alecja/electrochemistry_sim . Looking forward, we aim to apply this methodologyto study the various parameter regimes of electrochemical voltammetry experiments where13 b) (a) FIG. 2. Example CVs using the above method for at various scan rates, for the case with (a)no adsorption and (b) with adsorption. The units for current here are
InF A , while V is in volts.Mathematically, “no adsorption” means simulations using only Eq. 9 and not Eq. 10, i.e. ν ads and ν des are set to zero. “With adsorption” implies using both Eqs. 9 and 10, with k solf and k solb setto zero. The parameters for these simulations are D = 5 . · − cm /sec, T = 300 K, N T = 500, k = 10 cm/sec, α = 0 . adsorption plays a key role. After all, the addition of another timescale/parameter set cor-responding to adsorption/desorption rates, plus the potential asymmetry in redox behaviorof charge/uncharged species on a surface, implies that a realistic description of electrochem-ical dynamics and kinetics can and must go beyond the standard Nernstian and irreversibleregimes that we have all grown up with. 14 V. ACKNOWLEDGEMENTS
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, and the National Science Foun-dation under Grant No. DMS-2012286.
V. DATA AVAILABILITY
The data that support the findings of this study are available within the article.
VI. APPENDIX I: DERIVATION OF EQUATION 17
Below we will carry out a general derivation of Eq. 17 from Eq. 16, for either species A or B , with arbitrary initial conditions. We first break the integrand into two pieces, f ( s ) − f ( t ) √ t − s and f ( t ) √ t − s , where f ( s ) = b ( s,
0) + a ( s ) c j ( s, c j ( t,
0) = c j (0 ,
0) + ( D j π ) − . [ (cid:90) t [ b ( s,
0) + a ( s ) c j ( s, − [ b ( t,
0) + a ( t ) c j ( t, √ t − s ds + (cid:90) t [ b ( t,
0) + a ( t ) c j ( t, √ t − s ds ] c j ( t,
0) = c j (0 ,
0) + ( D j π ) − . [∆ s ( N − (cid:88) i = i b ( s i ,
0) + a ( s i ) c j ( s i , − b ( t, − a ( t ) c j ( t, √ t − s i + b (0 , − b ( t,
0) + a (0) c j (0 , − a ( t ) c j ( t, √ t ) + 2 (cid:112) ( t )( b ( t,
0) + a ( t, c j ( t, . (18)Rearrangement of all c j ( t,
0) terms to the left hand side yields c j ( t, (cid:112) D j π − a ( t ) √ t + ∆ s ( N − (cid:88) i =1 a ( t ) √ t − s i + a ( t )2 √ t )] − c j (0 , (cid:112) D j π =∆ s ( N − (cid:88) i =1 b ( s i ,
0) + a ( s i ) c j ( s i , − b ( t, √ t − s i + 0 . b (0 , − b ( t,
0) + a (0) c j (0 , √ t ) + 2 b ( t, √ t. (19)15inally, isolation of c j ( t,
0) yields the final result, c j ( t,
0) = ∆ s ( N − (cid:80) i =1 b ( s i , a ( s i ) c j ( s i , − b ( t, √ t − s i + b (0 , − b ( t, a (0) c j (0 , √ t ) + 2 b ( t, √ t + c j (0 , (cid:112) D j π (cid:112) D j π − a ( t )[2 √ t − ∆ s ( N − (cid:80) i =1 1 √ t − s i + √ t )] (20)Note that Eq. 17 is recovered exactly when the index j is replaced with B and the initialcondition c B (0 ,
0) = 0 is substituted.
REFERENCES N. Elgrishi, K. J. Rountree, B. D. McCarthy, E. S. Rountree, T. T. Eisenhart, and J. L.Dempsey, J. Chem. Educ. , 197 (2018). M. Bogdan, D. Brugger, W. Rosenstiel, and B. Speiser, J. Cheminform. (2014). E. P. Randviir, Electrochim. Acta , 179 (2018). A. J. Bard and L. R. Faulkner,
Electrochemical Methods: Fundamentals and Applications (Wiley, New York, 2001). M. Van Soestbergen, Russ. J. Electrochem. , 570 (2012). A. M. Limaye and A. P. Willard, J. Phys. Chem. C , 1352 (2020). C. Lin, X. Jiao, K. Tschulik, C. Batchelor-Mcauley, and R. G. Compton, J. Phys. Chem.C , 16121 (2015). A. J. Samin, AIP Adv. (2016). M. Yang and R. G. Compton, J. Electroanal. Chem. , 68 (2019). M. Yang and R. G. Compton, J. Phys. Chem. C , 18031 (2020). T. V. Shea and A. J. Bard, Anal. Chem. , 2101 (1987). A. W. Bott, Curr. Sep. , 9 (1999). H. Chen, J. R. Elliott, H. Le, M. Yang, and R. G. Compton, J. Electroanal. Chem. (2020). R. Guidelli, R. G. Compton, J. M. Feliu, E. Gileadi, J. Lipkowski, W. Schmickler, andS. Trasatti, in
Pure Appl. Chem. (2014). W. Schmickler and E. Santos,
Interfacial Electrochem. (2010). I. Langmuir, J. Am. Chem. Soc. , 1361 (1918). A. N. Frumkin, J. Electroanal. Chem. , 152 (1964).16 A. J. Coffman, A. K. Harshan, S. Hammes-Schiffer, and J. E. Subotnik, J. Phys. Chem.C , 13304 (2019). A. C. Hindmarsh, P. N. Brown, K. E. Grant, S. L. Lee, R. Serban, D. E. Shumaker, andC. S. Woodward, ACM Trans. Math. Softw. , 363 (2005). H. Risken,
The Fokker-Planck Equations, Method of Solution and Applications , 2nd ed.(Springer-Verlag, 1989) p. 103. A. J. Coffman, W. Dou, S. Hammes-Schiffer, and J. E. Subotnik, J. Chem. Phys. ,234108 (2020). F. John,