Electrochemical kinetics at electrode–electrolyte interfaces limit the performance of devices including fuel cells and batteries. While the importance of moving beyond Butler–Volmer kinetics and incorporating the effect of electronic density of states of the electrode has been recognized, a unified framework that incorporates these aspects directly into electrochemical performance models is still lacking. In this work, we explicitly account for the density functional theory-calculated density of states numerically in calculating electrochemical reaction rates for a variety of electrode–electrolyte interfaces. We first show the utility of this for two cases related to Li metal electrodeposition and stripping on a Li surface and a Cu surface (anode-free configuration). The deviation in reaction rates is minor for cases with flat densities of states such as Li, but is significant for Cu due to nondispersive *d*-bands creating large variation. Finally, we consider a semiconducting case of a solid-electrolyte interphase consisting of LiF and Li_{2}CO_{3} and note the importance of the Fermi level at the interface pinned by the redox reaction occurring there. We identify the asymmetry in reaction rates as a function of discharge/charge naturally within this approach.

## I. INTRODUCTION

Reactions at electrode–electrolyte interfaces control the limits of operation of various electrochemical devices. In the case of batteries, interfacial kinetics control the rate of reactions at the anode and cathode, which ultimately leads to the limits on fast discharge and charge capability of modern Li-ion batteries.^{1,2} Fast charging is an important requirement for electric vehicles,^{3} while fast discharge is necessary for electric vertical take-off and landing aircraft.^{4} In particular, there is significant interest in modifying the reaction rates of electrodeposition and stripping of Li metal electrodes by modifying the nature of the so-called solid electrolyte interphase (SEI).^{5}

The typical kinetic rate law used for Li metal electrodes is Butler–Volmer.^{6} Recently, Boyle *et al.*^{7} highlighted the importance of moving beyond Butler–Volmer kinetics and incorporating the effects of reorganization within the Marcus–Hush–Chidsey (MHC) formalism^{8,9} within the low overpotential limit. In an accompanying paper for this special issue, we extend that analysis and show the importance of accounting for the full MHC kinetic rate behavior.^{10} However, most treatments assume that the electrode density of states (DOS) is constant and independent of the energy states/overpotential. In the conventional treatment, this approach has no explanatory power for any asymmetry observed in the Tafel response. Asymmetric Marcus Theory^{11,12} has sought to address this issue by arguing that force constants (that is, the curvature of Gibbs free energy vs reaction coordinate) differ between reactants and products. However, as we show, in many cases, the impact of the DOS’s energy dependence on predicted reaction rates can also be large.

In this work, we incorporate the effect of the DOS at the electrode using density functional theory (DFT) calculations directly and numerically show its importance in a variety of cases. The approach can be thought of as a generalization of Gerischer’s model^{13} expressed in the “language” of MHC. We begin by discussing the effect of the DOS for Li electrodeposition and stripping for Li metal electrodes. Given the relatively flat nature of the DOS of Li metal, it affects the kinetic rates to a lesser extent. Next, we discuss the case of metal electrodeposition on Cu. Given the DOS associated with the more localized *d*-band states, there is significant deviation between the kinetic behavior when accounting for its variation compared to assuming a constant DOS. This highlights that the electronic state of the metal is an important factor in determining the high-rate performance of anode-free cells.^{14} Finally, we discuss the case of an electrochemical redox reaction with the density of states determined by SEI components, for example, LiF and Li_{2}CO_{3}.

We also provide both our full analysis code open-source on GitHub as well as a user-friendly interface for visualization of the effect of different model parameters on the results presented herein that readers can view online.

## II. METHODS

### A. Theoretical approach

The insight of Chidsey^{9} that it is critical to consider the occupation of electronic states when assessing rate constants within a Marcus-type theory was important. Building on this insight, the rate expressions developed in MHC theory (we here adapt notation from Ref. 15) are

for the reductive and oxidative directions, respectively, where *V*($\epsilon \u2009$) denotes a coupling constant between the (presumed localized) states on the electrolyte molecule and the (presumed plane-wave-like) states in the solid electrode. This coupling constant includes the density of states as well as an overlap integral term. *λ* is the reorganization energy, *e* is the electronic charge, *η* is the applied overpotential, *k*_{B} is Boltzmann’s constant, and *T* is the absolute temperature. All energies are measured relative to the Fermi level *E*_{f}. *f*_{FD}, denoting the Fermi-Dirac distribution for the occupation of state by electrons, is given by

and its complement, (1 − *f*_{FD}), is the distribution for holes.

In prior analyses, either explicitly or implicitly, the |*V*($\epsilon \u2009$)|^{2} term is presumed to have weak energy dependence and brought outside the integral. This allows the simplification presented in Ref. 16, taking advantage of the fact that the integral is over the whole real line, and the Fermi–Dirac distributions for electrons and holes can be equivalently thought of as complements of each other, or equal but with the energy coordinate running in the opposite direction,

However, this treatment neglects the electrode’s DOS, which generally *does* have strong variations with energy, especially over scales larger than a few *k*_{B}*T*. The most extreme example of this is in a semiconductor, where the value is identically zero at and near the Fermi level. This, as we will show, can lead to dramatically different rate constant predictions. However, even in the case of a transition metal (which often has large spikes in the DOS due to highly localized *d*-bands), this assumption can break down.

In our analysis, we consider adapted versions of Eqs. (1) and (2), where we explicitly account for the electrode DOS $D(\epsilon \u2009)$,

Because DOSs are not defined over an infinite range of energies, we must truncate the integrals in our calculations. In practice, this makes no discernible difference in the predicted rate constants because the Gaussian terms fall off dramatically within a relatively small energy range from *E*_{f}.

We also note that the dependence of the density of states (and to a lesser extent, the coupling constant that is also part of the |*V*($\epsilon \u2009$)|^{2} term in Eqs. (1) and (2)) has been described before,^{17} but is generally either assumed to be weak, or is modeled with simple functional forms that presume, e.g., spherical Fermi surfaces or other approximations that may or may not be valid across the wide diversity of relevant surfaces. To our knowledge, this is the first analysis that explicitly accounts for the DFT-calculated DOS numerically in directly calculating rates for electrochemical reactions at electrode–electrolyte interfaces.

### B. Computational approach

In this work, we use the Julia programming language to calculate the integrals described above using the adaptive Gauss–Kronrod quadrature^{18} as implemented in the QuadGK.jl package.

Density functional theory calculations for Li and Cu were conducted using the GPAW package^{19} through the Atomic Simulation Environment package.^{20} Ion–electron interactions were treated using the projector augmented wave approach.^{21} For all calculations, a grid spacing of 0.16 Å and a 4 × 4 × 1 Monkhorst–Pack k-mesh were used.^{22} To improve self-consistent field convergence, Fermi smearing was applied to electron occupation with a width of 0.05 eV. All relaxations and analysis, unless otherwise specified, were conducted using the Bayesian error estimation functional with van der Waals correction (BEEF-vdW) exchange-correlation functional. DOSs were calculated with a 100 meV Gaussian smearing width.

Calculations for LiF/Li_{2}CO_{3} were performed in Quantum Espresso according to methods described in Ref. 23.

## III. RESULTS AND DISCUSSION

### A. Li metal anode

The first case is the redox processes occurring in a Li metal electrode. The overall redox reaction is given by

where the forward reaction denotes electrodeposition (charging), and the backward reaction denotes stripping (discharging). The reaction rates are dependent on the organic electrolyte used,^{7} and the analysis is presented in the case of an EC:DEC electrolyte. Figure 1 shows experimental data and best-fit results for each model (results modeled with Eq. (4) are indicated in figure legends by “MHC” and from Eqs. (5) and (6) by “MHC + DOS”). Note that over the experimental overpotential range of ±0.2 V, the two models fit equally well. This is due to the fact that the DOS is relatively flat over this range. Over the full range of the plot, one can start to see variations between the models, but only modest ones.

This case illustrates why the importance of this more generalized analysis could have been missed previously—at relatively modest overpotentials, many DOSs are “close enough” to flat that predictions of the traditional MHC model and the one we present here differ little. We will now turn to cases where the two models show larger discrepancies.

### B. Cu current collector

We next examine the case of depositing Li on an anode-free configuration on top of a Cu current collector.^{14} This is given by

where * represents a site on a Cu surface. In this case, the electronic states of the Cu surface become important. The DOS (a) and predicted rate constants (b) for the Cu(111) surface are shown in Fig. 2. We chose Cu(111) as it has the most desirable nucleation characteristics among the low Miller index surfaces,^{14} in addition to the lowest surface energy.^{24} Note that due to the high density of states in the Cu *d*-bands starting ∼1.5 eV below the equilibrium *E*_{f}, a jump of roughly an order of magnitude in the rate constant is observed at this overpotential. An overpotential of this magnitude is entirely feasible, particularly in the context of the fast charge/discharge rates currently targeted for many applications.

It is worth noting that the scale of variation in DOS (i.e., an order of magnitude) is approximately equal to the change in the predicted rate constant. Indeed, one should expect this provided *λ* is much smaller than the overpotential window of interest, since *λ* essentially parametrizes a sliding Gaussian filter through which the DOS impacts *k* via the denominator of the first exponential terms in Eqs. (5) and (6).

### C. Effect of solid-electrolyte-interphase (SEI)

We will now consider the case where the Li electrode is covered with an inorganic solid electrolyte interphase due to the spontaneous reaction with the organic electrolyte typically used in Li-metal batteries. While the exact nature of the SEI remains an active area of research, it is widely acknowledged that some of the inorganic phases present include Li_{2}CO_{3}, LiF, Li_{2}O, etc.^{25} We consider a model SEI consisting of an interface between LiF and Li_{2}CO_{3}, building on the work of Pan *et al.*^{26} Here, we consider a redox reaction, e.g., Li deposition, occurring at the SEI component. Most SEI components are semiconducting or insulating by design and, hence, the kinetics for this semiconductor–electrolyte interface are quite different.

Figure 3 shows the results for this system. Because the DOS is identically zero throughout the ∼4 eV bandgap and the equilibrium *E*_{f} position of the pristine interface is at midgap, there is a discrepancy of as much as eight orders of magnitude between the predictions of the two models, even at low overpotentials, indicating the criticality of considering the energy dependency of the DOS in this context. This further substantiates the remark above regarding the scale of variation of the DOS—because it goes to zero (over a window of many *λ*) in the bandgap, this scale tends to infinity, and the scale of the discrepancy in *k* likewise diverges.

### D. Effects of model parameters and assumptions

Several parameters in the form of energy scales control the behavior of these models. Some of their behavior has already been alluded to above, and we examine these effects more closely here. First is the position of the equilibrium Fermi level *E*_{f}. Thus far, we have assumed it to be at the position predicted by DFT for a pristine surface. However, many system variables can influence this energy. Bulk or surface defects can shift *E*_{f}. In an electrochemical system, the redox reaction sets the Fermi level at the electrode–electrolyte interface and is an independent variable.^{27}

Figure 4(a) shows the effect of shifting *E*_{f} on the SEI case. In particular, we presume that an alternative redox couple (or equivalently, a particular voltage on the Li/Li^{+} scale) is chosen such that the equilibrium *E*_{f} is near the valence band maximum of the SEI. In this case, the behavior at low overpotentials exhibits a marked asymmetry: at negative *η* (pushing further into the valence bands), *k* rises sharply, while at positive *η* (pushing back into the bandgap), it drops exponentially until reaching midgap, whereupon it rises again as the conduction band minimum begins to have an impact. Such an asymmetry depending on the redox potential has been recognized^{28} and observed^{29} for Li–O_{2} batteries where discharge product, Li_{2}O_{2}, is a wide-band gap insulator similar to the LiF/Li_{2}CO_{3} interface considered here.

Figure 4(b) shows the effect of varying the reorganization energy *λ*. As mentioned previously, the mathematical effect of *λ* is primarily to act as a Gaussian smoothing filter on the density of states. However, an increased *λ* also increases the offset between the filters in the positive and negative directions, with the overall result being a slowly varying and slightly outward-shifted *k* vs *η* curve.

To better build the intuition for the impact of these parameters, we have built an interactive visualization, which can be viewed at https://nbviewer.jupyter.org/github/aced-differentiate/MHC_DOS/blob/master/dataviz.ipynb. It features the DOSs for all materials considered in this work and allows the viewer to modify *E*_{f}, *λ*, and the temperature and see in real time how rate constant predictions change.

It is worth noting that we have assumed that the electronic states of the electrode do not change along the reaction coordinate. This is exact for outer-sphere electron transfer reactions and is only approximate for inner-sphere electron transfer.^{30} For inner-sphere electron transfer, incorporating the electronic density of states along the reaction coordinate may become important. For the electrodeposition cases considered in this work, the approximation made is expected to be valid for the case of depositing Li on Li, while for the case of depositing Li on Cu may require re-visiting.

## IV. CONCLUSIONS

In this work, we explore the role of the electronic density of states at the electrode using density functional theory calculations and incorporate that numerically into a generalized Marcus–Hush–Chidsey kinetics formalism. We discuss three cases, electrodeposition on (i) Li surface and (ii) Cu surface, and (iii) redox reactions at an SEI component chosen here to be the LiF/Li_{2}CO_{3} interface. In the case of the Li surface, due to a relatively flat density of states, we find minor deviations from the standard MHC rate law. In the case of a Cu surface, mimicking an anode-free Li metal battery, the *d*-bands of Cu can significantly modify the rate constants, especially at high rates of charge and discharge. Finally, we show the drastic changes in reaction rates for a semiconducting interface and note the importance of the location of the Fermi level, which is pinned by the redox reactions occurring at that interface. We show strong asymmetric behavior in the SEI case, similar to that noted for Li_{2}O_{2} in Li–O_{2} batteries.

We believe this approach has broad promise for improved understanding of electrokinetics at solid–liquid interfaces such as those investigated here and could also be extended in the future to solid–solid interfaces such as those discussed in Refs. 31 and 32 by the inclusion of computed DOSs for the materials on both sides of the interface.

We have also made all analysis codes available open-source so that the community can take advantage of these models.

## ACKNOWLEDGMENTS

R.K. and V.V. gratefully acknowledge support from the Advanced Research Projects Agency Energy (ARPA-E) under Grant No. DE-AR0001211. This research was also supported by the Carnegie Mellon University Manufacturing Futures Initiative, which was made possible by the Richard King Mellon Foundation. We thank Yumin Zhang, Shang Zhu, and Zeeshan Ahmad for sharing structures used in this work.

## DATA AVAILABILITY

The data and code that support the findings of this study are openly available on GitHub at https://github.com/aced-differentiate/MHC_DOS.

## REFERENCES

_{ET}and i

_{STM}vs potential curves

_{2}O

_{2}in Li–O

_{2}batteries

_{2}O

_{2}