Electrodeposition and stripping are fundamental electrochemical processes for metals and have gained importance in rechargeable Li-ion batteries due to lithium metal electrodes. The electrode kinetics associated with lithium metal electrodeposition and stripping is crucial in determining the performance at fast discharge and charge, which is important for electric vertical takeoff and landing (eVTOL) aircraft and electric vehicles (EV). In this work, we show the use of Marcus–Hush–Chidsey (MHC) kinetics to accurately predict the Tafel curve data from the work of Boyle et al. [ACS Energy Lett. 5(3), 701 (2020)]. We discuss the differences in predictions of reorganization energies from the Marcus–Hush and the MHC models for lithium metal electrodes in four solvents. The MHC kinetic model is implemented and open-sourced within Cantera. Using the reaction kinetic model in a pseudo-2D battery model with a lithium anode paired with a LiFePO4 cathode, we show the importance of accounting for the MHC kinetics and compare it to the use of Butler–Volmer and Marcus–Hush kinetic models. We find significant deviation in the limiting currents associated with reaction kinetics for the three different rate laws for conditions of fast charge and discharge relevant for eVTOL and EV, respectively.

There is a need for improved energy density and specific energy for electrifying transportation and aviation.1,2 Lithium metal electrodes are a promising avenue to improve the energy density to meet these requirements.3 Reactions and their kinetics at electrode–electrolyte interfaces limit the performance of lithium metal electrodes.4 For a lithium metal electrode, during charge, lithium ions from the electrolyte (solid or liquid) deposit at the metal electrode, popularly known as electrodeposition, while during discharge, lithium metal is oxidized into lithium ions through a process known as stripping.5 

Fast charging is an important requirement for electric vehicles (EV),3 while fast discharge is an important requirement for the takeoff and landing segment of electric vertical takeoff and landing (eVTOL) aircraft.6 The morphological instabilities related to both fast discharge and charge have been well-documented leading to pitting and dendrite formation,7,8 thereby leading to poor cycle life.

Despite the enormous attention paid to understanding morphology, the kinetic rate behavior at the lithium metal electrode–electrolyte interface has received much less attention. Typical lithium metal models usually invoke Butler–Volmer (BV) kinetics,9 despite wide-recognition from the electrodeposition community that kinetic laws beyond Butler–Volmer are needed to describe fast discharge and charge conditions. Recently, Boyle et al.10 measured lithium electrodeposition and stripping kinetics by using transient voltammetry. Based on the measured data, they observed significant deviations from Butler–Volmer kinetics and argued for the use of Marcus–Hush (M–H) kinetics implemented as a low overpotential approximation of the Marcus–Hush–Chidsey (MHC) formalism. However, the low overpotential approximation breaks down for fast charge and discharge regimes.

In this work, building on the work of Boyle et al.,10 we show that the Marcus–Hush–Chidsey (MHC) model implemented as a uniformly valid closed form approximation developed by Zeng et al.11 agrees with the lithium electrodeposition and stripping data.10 We incorporate the MHC kinetic model into a pseudo-2D battery model and demonstrate for constant current fast charge and discharge applications up to a current density of 30 mA cm−2 and fast discharge eVTOL missions with discharge current densities of up to 26 mA cm−2. We find that the MHC kinetic formula is important to predict performance under high rate operating conditions in terms of limiting current densities and the distribution of charge transfer current through electrode thicknesses at fast discharge. To enable wide-spread use, we open-source the kinetic model for lithium metal into Cantera.12 We believe that future studies with lithium metal electrodes should utilize this formalism.

There are several approaches to modeling electron transfer kinetics,13 beginning in the 1930s with the Butler–Volmer model,14,15 microscopic electron transfer kinetics proposed by Marcus16 in the 1950s, kinetics of electrode processes examined by Hush,17,18 and, more recently, heterogeneous electron transfer explained using the Marcus–Hush–Chidsey19 kinetic model. Each of these kinetic models is presented with required modifications for their suitable interpretation in the context of electrochemical performance models of Li-ion batteries. Each kinetic model is presented for the following reaction:

Li++eLi.
(1)

The charge transfer kinetics at the interface of an electrode in electrochemical energy devices is often modeled using Butler–Volmer kinetics.13,20 Over the past several decades, the Butler–Volmer (BV) formalism has been used to model the kinetics at lithium metal anodes of Li-ion batteries.9,20–24 The BV formalism can be written as

kredBV=k0BV.expαFRTη,
(2)
koxBV=k0BV.exp(1α)FRTη,
(3)

where kred and kox are the reaction rates for reduction and oxidation, respectively, k0BV is the rate constant, η is the applied overpotential, and α is transfer coefficient that represents the position of the transition state of the reaction.25 The transfer coefficient, α, is a constant value that is generally 0.5 for single electron processes such as the one shown in [Eq. (1)]. The BV model posits a linear relationship between ln(kredBV/k0BV) and the overpotential η with a slope of −αF/RT, which is independent of the potential. However, experiments by Savéant and Tessier26 and many other studies show considerable deviations from linearity.13 The constant transfer coefficient independent of the potential is hence considered only a qualitative metric from BV theory and not universally applicable.13,25

Around 65 years ago, Marcus developed a kinetic model for microscopic outer-sphere homogeneous electron transfer.16 The theory has been used to model kinetics in several chemical and biological systems.27 The rate expressions for the Marcus model can be written as

kred/oxM=k0M.exp(ΔG±λ)24λ.RT,
(4)

where ΔG is the free energy change on reduction and λ is the reorganization energy, which required reorganizing the nuclear configuration of the reactants and solvent to the product state.13 One of the important predictions from Marcus theory is dependence of the transfer coefficient on the overpotential. The relevant derivations can be found elsewhere.13,16,26,28 The potential dependence of α is expressed as

α=12+Fη4λ
(5)

which can be incorporated by Eqs. (2) and (3). However, classical BV theory treats α as a constant.13 The use of a potential dependent α was shown by Savéant and Tessier26 and has also been used to explain recent lithium electrodeposition and stripping data by Boyle et al.10 

The Marcus kinetic model is developed for homogeneous electron transfer reactions where the two species involved in the reaction are dissolved in a solution. However, for applications such as metal electrodes, the charge transfer reaction is heterogeneous,25,29 i.e., between a metal electrode and a species dissolved in solution. For heterogeneous electron transfer with metal electrodes, the reaction kinetics will depend on electrons from different energy levels in the metal electrode conduction band.28 The dependence of reaction rates on the “electronic energy levels in metal” was noted by Marcus in 1965,28 and the incorporation of Fermi–Dirac statistics and its validation was shown by Chidsey in 1991,19 with a gold electrode and an electroactive ferrocene group.11,25,29 The rate equation used by Chidsey is generally referred to as the Marcus–Hush–Chidsey model,

kred/oxMHC=A.     exp(xλFη)24λ.RT.dx1+exp(x/RT),
(6)

where A is a pre-exponential factor that accounts for the density of metallic states and the electronic coupling strength.11,25 Compton and coworkers explained that the assumption of A being independent of potential has generally proven to be reasonable with experiments.19,25,30 Furthermore, for lithium, the density of states is relatively flat based on density functional theory calculations; the validity of this assumption is tackled in an accompanying paper for this Special Issue.31 The activation energy in Eq. (4) is recast using the electron energy, described as x, defined relative to the Fermi level.11,25,29

The net current density from the reaction, j, can be written as j = j0·(kred − kox), where j0 is the exchange current density. The net current for the BV model can be obtained using Eqs. (2) and (3),

jBV=j0BV.(kredBVkoxBV).
(7)

For applying Marcus theory, as noted previously, Boyle et al.10 employed a potential dependent α13,16,26,28 as shown in Eq. (5) within a BV current expression, as shown in Eq. (7). This model is referred to as the Marcus–Hush (M–H) kinetic model. The net current of the M–H model jM–H would have its exchange current density j0MH and follow the functional form shown in Eq. (7) along with α described by Eq. (5).

A closed form expression for the net current using the MHC model is required to avoid the evaluation of the MHC integral in Eq. (6), which is cumbersome to implement in a battery model. Bazant and co-workers have developed a uniformly valid closed form approximation of the integral,11 where they obtain a net current density expression using dimensionless reorganization energy, λ* = λ/RT, and dimensionless overpotential η* = F · η/RT,

jMHCS.πλ*.tanhη*2.erfcλ*1+λ*+η*22λ*
(8)

where S is used as a scaling factor for calculating the current density and erfc is the complementary error function. The closed form approximation provides the exchange current density to be11 

j0MHCS.πλ*2.erfcλ*1+λ*2λ*,
(9)

which is a function of the reorganization energy alone. This approximation has been shown to be consistent with the numerical integration of [Eq. (6)], with the accuracy increasing with the overpotential.11 Given that the goal of this work is to examine the implications of predictions from different kinetic models for high power applications of lithium batteries, we find that [Eq. (8)] provides the required flexibility and simplicity of the BV model while preserving the high accuracy of the MHC kinetic model.

We use transient voltammetry data from Boyle et al.10 collected using ultramicroelectrodes for LiPF6 in four solvents: (i) propylene carbonate (PC), (ii) diethyl carbonate (DEC), (iii) 1:1 by volume ethylene carbonate: diethyl carbonate (EC:DEC), and (iv) EC:DEC with 10% fluoroethylene carbonate (EC:DEC w. 10% FEC) to fit λ′s, which we will refer to as λM–H for the Marcus–Hush model and λMHC for the MHC model, respectively. The M–H model is implemented using Eq. (5), and the MHC model is implemented using Eq. (8). We also fit values for the exchange current densities j0M and j0MHC for the two kinetic models. The curve fitting was performed using the non-linear least squares method in MATLAB r2019b.

While the fundamental charge-transfer kinetic models described above are of, no doubt, great theoretical interest, do they result in a meaningful performance difference at technologically relevant scales? As shown in the aforementioned works,11,25 the models diverge only when the overpotential η is significant, relative to the solvent reorganization energy λ. With the push toward extreme fast charge (XFC) rates to facilitate EV consumer adoption and fast discharge for eVTOL, the probability of realizing suitably large overpotentials increases and merits investigation.

To explore the impact of detailed charge transfer kinetics on battery performance, we employ here a physically based pseudo-2D (P2D) model based on that developed by Doyle et al.9 The model simulates the performance of a porous LiFePO4–carbon (LFP/C) cathode, using porous electrode theory, paired with a dense lithium metal anode. The P2D model processes, geometry, and microstructure are illustrated in Fig. 1. We implement the three charge-transfer kinetic models as part of the open-source chemical kinetics software package Cantera,12 which is used to manage the thermo-kinetic calculations in a generalized manner.

FIG. 1.

Schematic illustration of P2D model processes. The model depicts a dense Li anode paired with a porous LFP–carbon cathode. At each electrode–electrolyte interface, faradaic charge transfer (jFar) across charged double layers is calculated according to three separate kinetic mechanisms: Marcus–Hush, Marcus–Hush–Chidsey, and Butler–Volmer. The simulations are run galvanostatically, with a user-specified external current (jext) boundary condition. The anode model includes an effective solid electrolyte interphase (SEI) model to incorporate the influence of the SEI current jSEI on overpotentials, but neglects many of the finer SEI details. In the cathode active material, intercalated lithium diffuses between the particle surface and interior at a rate NLi.

FIG. 1.

Schematic illustration of P2D model processes. The model depicts a dense Li anode paired with a porous LFP–carbon cathode. At each electrode–electrolyte interface, faradaic charge transfer (jFar) across charged double layers is calculated according to three separate kinetic mechanisms: Marcus–Hush, Marcus–Hush–Chidsey, and Butler–Volmer. The simulations are run galvanostatically, with a user-specified external current (jext) boundary condition. The anode model includes an effective solid electrolyte interphase (SEI) model to incorporate the influence of the SEI current jSEI on overpotentials, but neglects many of the finer SEI details. In the cathode active material, intercalated lithium diffuses between the particle surface and interior at a rate NLi.

Close modal

Despite the lower gravimetric capacity and electrical conductivity of LiFePO4, it has several benefits, such as lower material cost (6.3 Wh US$−1) than other standard cathode materials32 while having a high theoretical capacity among lithium metal phosphates, which tend to be more stable and have a longer cycling and calendar life.32 Additionally, advances in LFP cathodes have improved the specific capacity and electrical conductivity of the material,33,34 which, combined with favorable charge transfer kinetics at the electrode/electrolyte interface35 and high cyclability,34 make LFP a competitive cathode material for fast charge and high power applications.36,37

1. State variables

The solution variables tracked throughout the charge/discharge of the cell are as follows:

In the LiFePO4 phase,

  • XLi,ca is the intercalation fraction in the cathode (−) and

  • Φca is the electric potential of the cathode (V).

In the electrolyte phase,

  • Xk,elyte is the mole fraction of species k in the electrolyte (kmolk kmol−1) and

  • Φelyte is the electric potential of the electrolyte (V).

In the anode SEI,

  • ΦSEI/elyte is the electric potential in the SEI at the electrolyte interface (V).

In the dense lithium anode phase,

  • Φan is the electric potential of the lithium (V).

2. Model equations

The governing equations are derived from physically based conservation equations of mass, chemical elements, and electrical charge. The equations, described below, represent a set of differential-algebraic equations, which are integrated temporally to simulate galvanostatic charge–discharge with a specified external current jext. While the P2D model equations are common and well-represented throughout the literature,9,38–40 we present here a brief overview of the governing equations.

a. Cathode composition.

In the cathode, assuming spherical particles, diffusion only in the radial direction, and a constant diffusion coefficient DLi,ca, the governing equation for the intercalated lithium mole fraction is derived from conservation of mass and elements,

XLi,cat=(DLi,caXLi,ca)+ṡLiAintεcaCca,
(10)

where Cca is the total molar concentration of the cathode phase (mol m−3), and the term ṡLi is a molar source term due to the charge transfer reaction at the LFP particle surfaces (mol m−2 s−1). Aint is the specific surface area of the electrode/electrolyte interface (mint2m3). εca is the volume fraction of the active cathode phase. At the particle surface, it is evaluated via one of the three charge transfer kinetic models described above, with the relationship between the molar production rate and the Faradaic current given by Faraday’s law,

ṡk,elyte=νk,LijFarzF,
(11)

where νk,Li is the stoichiometric coefficient for intercalated lithium in the cathode phase (+1, here) and z represents the moles of charge transferred per mole of reaction (z = 1 in this study). A similar relationship is used for the electrolyte species production rates, below. Finally, ṡLiAint is equal to zero, internal to the particle.

b. Electrolyte composition.

Conservation of mass and elements is also applied to derive the governing equation for the electrolyte species mole fractions,

Xk,elytet=1CelyteεelyteNk,elyte+ṡk,elyteAint,
(12)

where Celyte is the electrolyte molar concentration (mol m−3, assumed to be constant for the incompressible electrolyte), εelyte is the electrolyte volume fraction, Nk, elyte is the molar flux (molk m−2 s−1), calculated according to dilute solution theory,41 and ṡk,elyte is the net production rate of electrolyte species k (molkmint2s1) due to heterogeneous reactions at the electrode/electrolyte interface. Within the electrolyte separator, the term ṡk,elyteAint equals zero.

c. Electric potentials.

Finally, governing equations for the electrode and electrolyte electric potentials are derived by applying charge conservation and assuming charge neutrality in bulk phase interiors. This yields an equation for the double-layer potential at all electrode–electrolyte interfaces ΔΦdl = Φelyte − Φed [V, where “ed” = “ca” (cathode) or “an” (anode)], and an algebraic constraint on the sum of all currents into the volume

ΦelyteΦedt=jdlCdl, edAint
(13)

and

0=jio+jel.
(14)

In these equations, jdl is the double-layer current per unit volume at the electrode/electrolyte interface (A m−3), Cdl,ed is the double-layer capacitance (F m−2), and jio and jel are the ionic and electronic current densities, respectively (A m−2, both defined as positive for the positive charge moving in the positive y direction). The double layer current is calculated by enforcing charge neutrality within the bulk of either the electrode or electrolyte phase (without loss of generality, we choose the electrode),

0=jel+jFar+jdl,
(15)

where jFar is the Faradaic current at the electrode/electrolyte interface due to charge transfer reactions, per unit volume (A m−3). The interfacial currents jFar and jdl are both defined such that positive j represents positive charge entering the electrode and leaving the electrolyte. Therefore, positive jdl represents an increase in the positive charge in the electrolyte side of the charged double layer.

The handling of the double layer here represents one potential point of difference, relative to other P2D models. Some models38,39 assume that the double layer responds infinitely quickly such that jdl is zero at all times. This converts the problem to an additional algebraic constraint such that the Φelyte must always be set to that value, which results in jFar = −∇jel. The approach in this study results in a differential equation and captures the dynamic response of the double layer during battery charge and discharge, which can have non-negligible impacts, particularly for XFC applications.

Finally, although Boyle et al. demonstrated that SEI dynamics do not significantly impact the electrode kinetics in their transient voltammetry experiments,10 SEI will certainly form and influence performance during charge–discharge. We estimate here the impact of an SEI with constant resistance (0.698 mΩ m2 is used here42). This imposes an additional algebraic constraint in that the current across the SEI (modeled here as Ohmic in nature) must equal the Faradaic current at the electrolyte/SEI interface. In other words, the model determines the SEI electric potential at the electrolyte interface (Φelyte/SEI) such that jFar = jSEI. The charge transfer kinetics use the potential difference Φelyte − ΦSEI/elyte to calculate the overpotential η = ΔΦ − ΔΦeq. The charge-transfer reaction here takes place at the electrolyte-SEI interface, where the solvent molecules must re-arrange to allow for transfer of Li+ into the SEI. Charge-transfer at the SEI/anode interface is assumed to have negligible resistance. Clearly, the actual electrochemistry and transport in the SEI are much more complex. As demonstrated below, charge-transfer at the Li anode is not rate-limiting at the conditions explored herein, and the present first-order approximation is sufficient for present purposes (demonstrating the impact of varying kinetic mechanisms on P2D results).

d. Model and simulation parameters.

Model parameters were taken from previous experimental and modeling studies for Li–LiFePO4/C.10,33,41–44 The P2D simulations are presented here to demonstrate the impact of charge transfer kinetics within a relevant application, not as a fully validated, predictive modeling study. As such, model parameters were adapted to qualitatively match rate capability data from the recent literature45 but have not been rigorously fit to the data. This is left as an extension for a future study. Validation results and a detailed list of model parameters are presented in the supplementary material.

Fundamental electrochemical data and molecular theory support MHC kinetics as most appropriate for both the LFP/C cathode and Li metal anode. The data are not readily fit by BV kinetics, and the presence of electron conductors (also the lack of any experimentally observed “inverted region”) suggest that Marcus–Hush kinetics are not the most appropriate.10,43 Regardless, it is worth exploring the comparative performance of the three mechanisms, in the context of the P2D model, to understand the magnitude of the errors made by researchers adopting either the commonly employed M–H or BV kinetics for Li anodes of LFP/C cathodes in device simulations.

Simulations were run for two conditions: (i) galvanostatic charge–discharge at current densities ranging from 1 mA cm−2 to 30 mA cm−2 for cathode thicknesses ranging from 20 μm to 60 μm and (ii) simulated eVTOL flight powered by a battery with a 40 μm thick cathode, including 1-min takeoff and landing segments at 25 mA cm−2 and 26 mA cm−2, respectively, and a 2-min cruise at 10 mA cm−2.

The source code for the model, including instructions for installation and running, and all parameter values are made available via GitHub.46 Although we plan to incorporate the charge transfer kinetics into an official future release of Cantera, currently, the required functionality can be obtained by using a forked version of the chemical thermo-kinetics software, published as a GitHub repository.47 

The transient voltammetry data for lithium electrodeposition and stripping from Boyle et al.10 is evaluated for the M–H model and the MHC model described in Sec. II A to find the exchange current densities, j0’s, and the reorganization energies, λ′s. We compare the Marcus–Hush (M–H) kinetic model proposed by Boyle et al. and our implementation of the MHC model along with the BV model. Each of these is evaluated in the context of using these kinetic models within lithium battery models to supplant the existing BV models. We find that the MHC model is more appropriate to use within electrochemical battery models. We compare the predictions from the P2D model in fast charge and discharge conditions for the BV, M–H, and MHC kinetic models.

The best fit exchange current densities for the M–H model (j0MH) and the MHC model (j0MHC), along with the reported linear fit values (j0L.fit) from Boyle et al.,10 are shown in Table I. The values of exchange current density fit for the M–H and MHC model in this work agree closely with each other. However, the j0L.fit values are seen to be marginally higher than j0MH and j0MHC. The trends in exchange current density hold in each of the estimates, where the exchange current density of PC < DEC < EC:DEC < EC:DEC wt. 10% FEC. As noted by Boyle et al., the current density predicted that the transient voltammetry data are not affected by the solid electrolyte interphase (SEI),10 and the reaction is in kinetic control.

TABLE I.

The exchange current density values as obtained from the M–H model, joMH, and the MHC model, joMHC, and values reported by Boyle et al. using a linear fit.

Solvent j o M H ( m A c m 2 ) joL.fit(mAcm2) Ref. 10  j o M H C ( m A c m 2 )
PC  1.9  2.6  1.9 
DEC  2.2  3.7  2.2 
EC:DEC  8.8  10.4  8.6 
EC:DEC wt. 10% FEC  14.5  16.0  13.8 
Solvent j o M H ( m A c m 2 ) joL.fit(mAcm2) Ref. 10  j o M H C ( m A c m 2 )
PC  1.9  2.6  1.9 
DEC  2.2  3.7  2.2 
EC:DEC  8.8  10.4  8.6 
EC:DEC wt. 10% FEC  14.5  16.0  13.8 

The best fit values for the reorganization energy for the M–H model (λM–H) and the MHC model (λMHC) along with the reported M–H fit values from Boyle et al.10 are shown in Table II. We observe that the λM–H values of our study and λM–H⋄ from Boyle et al. are similar; however, λMHC values are seen to be generally lower than λM–H values for each of the solvents. The relative values and the trend for the reorganization energy between the four solvents hold across M–H and MHC models. The calculation of λM–H and λM–H⋄ neglects the distribution of electrons, while λMHC calculated using Eq. (8) implicitly accounts for it. While the exchange current density values of our M–H and MHC models agree closely, the reorganization energy values are consistently lower for MHC by about 0.12 ± 0.1 eV.

TABLE II.

The reorganization energy values as obtained from the M–H model, λM–H, and the MHC model, λMHC, and values reported by Boyle et al. from another implementation of the M–H model.

Solvent λM–H(eV), [95% CI] λM–H⋄(eV) Ref. 10  λMHC(eV), [95% CI]
PC  0.33 [0.30,0.36]  0.34  0.21 [0.18,0.24] 
DEC  0.38 [0.29,0.46]  0.34  0.25 [0.18,0.24] 
EC:DEC  0.34 [0.31,0.37]  0.3  0.22 [0.19,0.26] 
EC:DEC wt. 10% FEC  0.31 [0.29,0.33]  0.28  0.19 [0.17,0.21] 
Solvent λM–H(eV), [95% CI] λM–H⋄(eV) Ref. 10  λMHC(eV), [95% CI]
PC  0.33 [0.30,0.36]  0.34  0.21 [0.18,0.24] 
DEC  0.38 [0.29,0.46]  0.34  0.25 [0.18,0.24] 
EC:DEC  0.34 [0.31,0.37]  0.3  0.22 [0.19,0.26] 
EC:DEC wt. 10% FEC  0.31 [0.29,0.33]  0.28  0.19 [0.17,0.21] 

The MHC model consistently predicts lower values of reorganization energy for all four solvents while consistently maintaining similar values for exchange current density with the M–H model. The only difference between the M–H model (implemented using a potential dependent transfer coefficient α) and the MHC model [implemented using the uniformly valid closed form approximation in Eq. (8)] is the use of Fermi–Dirac statistics to account for distribution of electrons in different energy levels. Hence, the difference in reorganization energies can be attributed to the inclusion of the distribution of electrons over different energy levels. Furthermore, the M–H implementation discussed here and as implemented in the work of Boyle et al. is a low overpotential approximation of the MHC model. Hence, it follows that the use of the M–H model should be limited to overpotential regions where η < λ. Furthermore, Bai and Bazant43 have confirmed from first principles calculations that the reorganization energies predicted by the MHC model are more accurate compared to typical Marcus or M–H models. Bai and Bazant43 also noted that Marcus or M–H models require larger reorganization energies to fit electrode kinetics data because of the inverted region predicted for η > λ. Between the two models, the MHC rate expression implemented using the closed form approximation in Eq. (8) is more appropriate to use within an electrochemical battery model, as shown in discussions pertaining to Figs. 2 and 3.

FIG. 2.

The predictions of current density as provided by the BV model, the M–H model (proposed by Boyle et al.),10 and the MHC model for four different solvents, namely, (a) PC, (b) DEC, (c) EC:DEC, and (d) EC:DEC wt. 10% FEC. The exchange current densities and the reorganization energy values for each of the fits are shown in Tables I and II, respectively. The root mean square error (RMSE) in mA/cm2 for M–H and MHC models is very similar, while the goodness of fit characterized by R2 is the same for both M–H and MHC for each solvent (R2 for both M–H and MHC models are 0.997, 0.987, 0.992, and 0.997 for PC, DEC, EC:DEC, and EC:DEC wt. 10% FEC, respectively). It is evident that both M–H and MHC models explain the data very well and collapse on each other at low overpotentials; however, at overpotentials greater than an absolute value of 0.25 V, the difference in predictions is significant. This is because the M–H model posits an inverted region for current at overpotentials η > λ, while the MHC model predicts a plateau. Given that these models are examined in the context of using them in lithium battery models under high power conditions, it is important to use the more appropriate model for overpotentials greater than 0.25 V. Predictions from the BV model are also shown here for reference, where the BV model uses the same exchange current density of the M–H model. As seen for each solvent, the BV model deviates significantly for overpotentials greater than 0.1 V.

FIG. 2.

The predictions of current density as provided by the BV model, the M–H model (proposed by Boyle et al.),10 and the MHC model for four different solvents, namely, (a) PC, (b) DEC, (c) EC:DEC, and (d) EC:DEC wt. 10% FEC. The exchange current densities and the reorganization energy values for each of the fits are shown in Tables I and II, respectively. The root mean square error (RMSE) in mA/cm2 for M–H and MHC models is very similar, while the goodness of fit characterized by R2 is the same for both M–H and MHC for each solvent (R2 for both M–H and MHC models are 0.997, 0.987, 0.992, and 0.997 for PC, DEC, EC:DEC, and EC:DEC wt. 10% FEC, respectively). It is evident that both M–H and MHC models explain the data very well and collapse on each other at low overpotentials; however, at overpotentials greater than an absolute value of 0.25 V, the difference in predictions is significant. This is because the M–H model posits an inverted region for current at overpotentials η > λ, while the MHC model predicts a plateau. Given that these models are examined in the context of using them in lithium battery models under high power conditions, it is important to use the more appropriate model for overpotentials greater than 0.25 V. Predictions from the BV model are also shown here for reference, where the BV model uses the same exchange current density of the M–H model. As seen for each solvent, the BV model deviates significantly for overpotentials greater than 0.1 V.

Close modal
FIG. 3.

A comparison of the predictions for current density from M–H and MHC models. (a) shows the difference between jMHC and jM–H for the potential window between −0.5 V and 0.5 V. (b) The difference between jMHC and jM–H for the potential range of −0.25 V and 0.25 V, which are the limits of the overpotential in the dataset. We observe that the difference is about 50 mA/cm2 for PC and DEC, while it goes up to 200 mA/cm2 for EC:DEC at 0.5 V of the overpotential and about 250 mA/cm2 at 0.5 V for EC:DEC wt. 10% FEC. Within an overpotential region of 0.25 V, the greatest error is about 1.5 mA/cm2 for EC:DEC wt. 10% FEC. when the absolute value of current density is over 200 mA/cm2, which is a difference of under 1%. As explained in Fig. 2, it is clear that both the M–H and MHC models explain the data with almost equal accuracy, with low error and small difference between the two models for overpotential regions within an absolute value of 0.25 V; however, the difference in predictions is significant when we consider a slightly higher overpotential window of about 0.5 V.

FIG. 3.

A comparison of the predictions for current density from M–H and MHC models. (a) shows the difference between jMHC and jM–H for the potential window between −0.5 V and 0.5 V. (b) The difference between jMHC and jM–H for the potential range of −0.25 V and 0.25 V, which are the limits of the overpotential in the dataset. We observe that the difference is about 50 mA/cm2 for PC and DEC, while it goes up to 200 mA/cm2 for EC:DEC at 0.5 V of the overpotential and about 250 mA/cm2 at 0.5 V for EC:DEC wt. 10% FEC. Within an overpotential region of 0.25 V, the greatest error is about 1.5 mA/cm2 for EC:DEC wt. 10% FEC. when the absolute value of current density is over 200 mA/cm2, which is a difference of under 1%. As explained in Fig. 2, it is clear that both the M–H and MHC models explain the data with almost equal accuracy, with low error and small difference between the two models for overpotential regions within an absolute value of 0.25 V; however, the difference in predictions is significant when we consider a slightly higher overpotential window of about 0.5 V.

Close modal

A comparison between the different models and how they explain the transient voltammetry data is shown in Fig. 2. One of the first observations that can be made is the significant deviation of the BV model from the data for all solvents. This observation has been noted for several metal electrodeposition and stripping studies48,49 and pointed out for the case of lithium as well.10 Furthermore, we can see that the M–H model proposed by Boyle et al. and the MHC model from this work explain the data with similar accuracy; however, the current density predictions from the two models begin to deviate significantly after about 0.3 V of the overpotential. The M–H model predicts an inversion of current similar to the inverted region proposed by Marcus theory.16,28 As noted previously, the M–H implementation shown here and in Boyle et al. is a low overpotential approximation of the MHC model, while the MHC model implemented in this work is a uniformly valid approximation. Hence, we see that the two models are indistinguishable at lower overpotentials. The difference between the models is significant at overpotentials greater than 0.25 V and deviation is larger for larger overpotentials.

The deviation of the M–H model from MHC model starts at about 0.25 V of the overpotential. Figures 3(a) and 3(b) show the difference between the absolute values of current from the two kinetic models. Figure 3(b) shows the deviation within an overpotential region of 0.25 V, which is the maximum overpotential in the transient voltammetry data.10 We observe that, for the overpotential region within 0.25 V, the maximum deviation of about 1.5 mA/cm2 is caused for the high exchange current density solvent EC:DEC wt. 10% FEC at 0.25 V of the overpotential. The difference between the currents from the two kinetics models is much lower for the other three solvents. These differences are very small compared to the absolute value of the current at this potential, which is about 250 mA/cm2. On the other hand, on examining the difference in currents for the overpotential region within 0.5 V, we observe that, for EC:DEC wt. 10% FEC, the deviation of the prediction of current from the M–H model is about 250 mA/cm2. Similarly, for lower exchange current density solvents such as PC and DEC, the deviation is about 50 mA/cm2. Due to the inverted region predicted by the M–H model, the current predictions are always an underestimate at overpotentials greater than 0.25 V for the M–H model. We also observe an asymmetry in the electrodeposition and stripping data, which is consistent across all the solvents. To explore the asymmetry further, we fit the data using the asymmetric Marcus–Hush model,50 and the results are discussed in the supplementary material.

Based on the discussions here and the results shown in Figs. 2 and 3, while the MHC model collapses onto the M–H model at low overpotentials, the deviation of the predictions from the M–H model are substantial at slightly higher overpotentials here. Since these deviations are significant and because the M–H model is only a low overpotential approximation of the MHC model, the MHC model implemented using the uniformly valid closed form approximation11 is better suited for use with electrochemical models of lithium batteries. The exchange current density fit for both models are identical, but the reorganization energies are lower by about 0.12 eV for MHC kinetic fits. Further experiments beyond overpotentials of 0.25 V are required to confirm the nature of the Tafel curves for lithium electrodeposition and stripping at higher overpotentials.

Implementing the kinetic models in Sec. II A within the P2D model framework from Sec. II B demonstrates the importance of accurate charge-transfer kinetics at the application scale. The predicted battery performance is controlled by a combination of electrochemical kinetic and transport phenomena in the electrolyte and solid cathode phases, each of which can play limiting roles under different conditions. The model results presented below identify conditions at which the Butler–Volmer and Marcus–Hush models predict qualitatively different battery behavior, compared to the appropriate Marcus–Hush–Chidsey kinetics.

As suggested by Fig. 2, the Li anode overpotentials predicted by the three kinetic mechanisms for any given current density deviate significantly for currents ≥40 mA cm−2, which are unlikely even for current XFC and other high-power applications. However, overpotentials for the three mechanisms deviate at much lower current densities at the LFP/C electrode,43 leading to important discrepancies in predicted battery performance. Specifically, for kinetically limited operating regimes (thin cathodes or large primary particle sizes), MHC kinetics predict limiting currents not captured by BV kinetic predictions. M–H kinetics, on the other hand, predict very low limiting currents, in contrast to the experimental results.41,45 Accurate model predictions will be critical to inform battery design and control strategies for cases such as XFC in EVs or eVTOL vehicles.

Figures 4 and 5 show predicted charge curves and available charge capacities for a Li–LFP/C cell operating at current densities ranging from 1 mA cm−2 to 30 mA cm−2 for cathode thicknesses of 20 μm, 40 μm, and 60 μm. As expected, the capacity decreases with increasing current and generally increases with increasing cathode thickness, with one exception. For the 60 μm thick cathode, we see a sharp drop in capacity due to electrolyte transport limitations in the cathode, which lead to a depletion of Li ions at the anode surface.

FIG. 4.

P2D model predictions for XFC with current densities from 1 mA cm−2 to 30 mA cm−2 for (a) 20 μm, (b) 40 μm, and (c) 60 μm cathodes. Results for the three kinetic models demonstrate the importance of accurate charge transfer kinetics at the cell level. For kinetically controlled conditions in the cathode (e.g., low cathode thickness and high current density), the MHC and M–H models predict lower limiting currents compared to the commonly used Butler–Volmer form. The M–H model is unable to sustain charging at current densities higher than roughly 1 mA cm−2–10 mA cm−2, due to the so-called “inverted region.” Note that a current density of 1 mA cm−2 corresponds to C rates of 1.2 C, 0.6 C, and 0.4 C for the 20 μm, 40 μm, and 60 μm thick cathodes, respectively.

FIG. 4.

P2D model predictions for XFC with current densities from 1 mA cm−2 to 30 mA cm−2 for (a) 20 μm, (b) 40 μm, and (c) 60 μm cathodes. Results for the three kinetic models demonstrate the importance of accurate charge transfer kinetics at the cell level. For kinetically controlled conditions in the cathode (e.g., low cathode thickness and high current density), the MHC and M–H models predict lower limiting currents compared to the commonly used Butler–Volmer form. The M–H model is unable to sustain charging at current densities higher than roughly 1 mA cm−2–10 mA cm−2, due to the so-called “inverted region.” Note that a current density of 1 mA cm−2 corresponds to C rates of 1.2 C, 0.6 C, and 0.4 C for the 20 μm, 40 μm, and 60 μm thick cathodes, respectively.

Close modal
FIG. 5.

P2D model predictions for XFC gravimetric capacity (mAh per g LFP) with current densities from 1 mA cm−2 to 30 mA cm−2 for (a) 20 μm, (b) 40 μm, and (c) 60 μm cathodes. Results demonstrate that the principle difference between the three mechanism capacities is the prediction of limiting currents. So long as i < ilim, the mechanisms predict nearly identical capacities. However, Marcus–Hush kinetics predict anomalously low limiting currents, while Butler–Volmer kinetics do not predict any limiting currents at all, both of which contradict experimental observations. For the 60 μm thick cathode at i = 30 mA cm−2, the low capacity is due to electrolyte transport limitations in the porous cathode, not due to kinetic limitations.

FIG. 5.

P2D model predictions for XFC gravimetric capacity (mAh per g LFP) with current densities from 1 mA cm−2 to 30 mA cm−2 for (a) 20 μm, (b) 40 μm, and (c) 60 μm cathodes. Results demonstrate that the principle difference between the three mechanism capacities is the prediction of limiting currents. So long as i < ilim, the mechanisms predict nearly identical capacities. However, Marcus–Hush kinetics predict anomalously low limiting currents, while Butler–Volmer kinetics do not predict any limiting currents at all, both of which contradict experimental observations. For the 60 μm thick cathode at i = 30 mA cm−2, the low capacity is due to electrolyte transport limitations in the porous cathode, not due to kinetic limitations.

Close modal

Comparing the three kinetic mechanisms, we see that the primary difference is in terms of the predicted limiting current. M–H kinetics predict rather low limiting currents at all thicknesses due to the inverted region. Specifically, M–H kinetics predict an inability to sustain current densities >1 mA cm−2 for the 20 μm thick cathode and >10 mA cm−2 for 40 μm and 60 μm cathodes, predictions that run counter to experimental observations.41,43–45 For currents where M–H kinetics predict battery operation, the charge profiles are nearly identical to those for Butler–Volmer and Marcus–Hush–Chidsey kinetics. However, clearly, the theory is not appropriate for the Li-LFP/C battery, given the availability of valence-band electrons in both electrodes.

The M–H model predictions in Figs. 4 and 5 are explained by the so-called “inverted region” in Fig. 2, where the current density decreases as a function of the overpotential η for η > λ. This phenomenon is shown in greater detail in Fig. 6, which plots the charging voltage for the M–H model for a 60 μm thick cathode charging at 20 mA cm−2. The voltage initially increases due to charging of the double layer. As the resulting Faradaic current density at the cathode–electrolyte interface increases commensurately with ΔΦdl, the difference between the jext and the total jFar decreases, and the slope of the voltage profile also decreases. However, once the overpotentials reach the point corresponding to jFar, max (roughly η = 0.35 V, in Fig. 2), further charging of the double layer decreases jFar, which leads to unconstrained growth of the double layer potential. This is observed in Fig. 6 at a capacity of roughly 0.1 μAh cm−2. As jFar goes to zero with further increases in ΔΦdl, the double layer current is essentially constant at jdl = jext, reflected in the linear growth for the cell potential at discharge capacities greater than 0.1 μAh cm−2 (note that capacity is proportional to time for galvanostatic discharge).

FIG. 6.

P2D-predicted charging profile for a 60 μm thick cathode with the M–H kinetic model and a charging current of 20 mA cm−2, demonstrating the transition from the charge transfer current predicted at low overpotentials to the “inverted region” at high overpotentials.

FIG. 6.

P2D-predicted charging profile for a 60 μm thick cathode with the M–H kinetic model and a charging current of 20 mA cm−2, demonstrating the transition from the charge transfer current predicted at low overpotentials to the “inverted region” at high overpotentials.

Close modal

Finally, it is worth noting that realistic initial conditions are crucial, when using M–H kinetics, due to the inverted region. Most kinetic models are “self-correcting” in that they can locate the equilibrium potential from any starting guess by simulating open circuit conditions. With M–H kinetics, however, if starting conditions place the kinetic mechanism in the inverted region, the double layer will continue to charge irreversibly, even at the zero-current open circuit condition, predicting an inability to charge or discharge the battery even at very low currents.

Although MHC kinetics do not predict an inverted region, they do predict limiting currents once barrierless charge transfer occurs. This is observed only for high current densities (30 mA cm−2) for the 20 μm thick cathode [Figs. 4(a) and 5(a)], where kinetics become limiting. In the 40 μm and 60 μm thick cathodes, simulations predict slight discrepancies between charge voltages using MHC and BV kinetics, on the order of 10 mV, for currents ≥20 mA cm−2. While MHC and BV kinetics therefore predict incredibly similar performance below the MHC limiting current, robust kinetic predictions over a range of operating conditions and cathode designs benefit from the use of MHC kinetics to correctly predict limiting currents. For the BV model, results for j ≥ jLim are qualitatively incorrect. Moreover, the different voltage profiles between MHC and BV kinetics result in different predictions of irreversible energy loss for the two mechanisms due to the discrepancy between the charging voltage and the equilibrium potential (i.e., the overpotential η). Integrating j × η as a function of time yields the irreversible energy loss during charging. Figure 7 shows the loss predicted by MHC kinetics, relative to the BV prediction. In other words, the figure represents the degree to which BV kinetics underestimate the irreversible energy loss during XFC on the order of 1 to 3 mWh gLFP1 at 20 to 30 mA cm−2.

FIG. 7.

Difference in predicted irreversible energy loss during battery charging, Marcus–Hush–Chidsey vs Butler–Volmer kinetic models. The energy loss is calculated for each mechanism by integrating the current density times the total overpotential. The plot shows the energy loss predicted by Marcus–Hush–Chidsey kinetics relative to that for the Butler–Volmer model. Because MHC predicts an inability to charge the 20 μm cathode at 30 mA cm−2, no entry is provided.

FIG. 7.

Difference in predicted irreversible energy loss during battery charging, Marcus–Hush–Chidsey vs Butler–Volmer kinetic models. The energy loss is calculated for each mechanism by integrating the current density times the total overpotential. The plot shows the energy loss predicted by Marcus–Hush–Chidsey kinetics relative to that for the Butler–Volmer model. Because MHC predicts an inability to charge the 20 μm cathode at 30 mA cm−2, no entry is provided.

Close modal

To explore the impact of kinetic models on fast discharge (relevant for high-power applications), Figs. 8 and 9 show discharge curves at 1 mA cm−2–30 mA cm−2 and a simulated eVTOL flight, respectively. The discharge curves in Fig. 8 are qualitatively similar to the corresponding charge curves in Fig. 4, in terms of both capacity trends and the comparison between mechanisms, but with overall lower capacities due to the slower LFP kinetics during discharge.43 As with charging, the primary difference between the mechanisms lies in the prediction of limiting current densities. For j < jLim, the mechanisms predict nearly identical results. This is also reflected in the simulated eVTOL flight in Fig. 9. Outside of the very beginning of takeoff and the very end of landing, the two mechanisms (MHC and BV) predict nearly identical potentials (therefore power densities). The subplots show the flight’s very beginning and ending in greater detail, where predicted potentials from the two mechanisms vary by as much as 40 mV during takeoff and 50 mV during landing, corresponding to power density discrepancies of 1.0 mW cm−2 and 1.3 mW cm−2, respectively.

FIG. 8.

P2D model predictions for battery discharge with current densities from 1 mA cm−2 to 30 mA cm−2 for (a) 20 μm, (b) 40 μm, and (c) 60 μm cathodes. For kinetically controlled conditions in the cathode (e.g., low cathode thickness and high current density), the MHC and M–H models predict limiting current densities, whereas the commonly used Butler–Volmer form does not. A current density of 1 mA cm−2 corresponds to C rates of 1.2 C, 0.6 C, and 0.4 C for the 20 μm, 40 μm, and 60 μm thick cathodes, respectively.

FIG. 8.

P2D model predictions for battery discharge with current densities from 1 mA cm−2 to 30 mA cm−2 for (a) 20 μm, (b) 40 μm, and (c) 60 μm cathodes. For kinetically controlled conditions in the cathode (e.g., low cathode thickness and high current density), the MHC and M–H models predict limiting current densities, whereas the commonly used Butler–Volmer form does not. A current density of 1 mA cm−2 corresponds to C rates of 1.2 C, 0.6 C, and 0.4 C for the 20 μm, 40 μm, and 60 μm thick cathodes, respectively.

Close modal
FIG. 9.

Predicted discharge curve of a time-dependent current simulating the changing load during takeoff, flight, and landing that an eVTOL might experience. The MHC and BV kinetics models predict nearly identical voltage and power throughout the flight time. The lone exception is early during takeoff and late during landing when the two predicted voltage profiles differ by roughly 50 mV (corresponding to a power density difference of roughly 1 mW cm−2).

FIG. 9.

Predicted discharge curve of a time-dependent current simulating the changing load during takeoff, flight, and landing that an eVTOL might experience. The MHC and BV kinetics models predict nearly identical voltage and power throughout the flight time. The lone exception is early during takeoff and late during landing when the two predicted voltage profiles differ by roughly 50 mV (corresponding to a power density difference of roughly 1 mW cm−2).

Close modal

For both charge and discharge, then, the three kinetic mechanisms predict nearly identical performance for all current densities below the limiting value for each mechanism. Even for currents just 1 mA cm−2 or 2 mA cm−2 below jLim (not shown here), predicted voltage discrepancies are roughly 10 mV. This is perhaps surprising, given differences of up to 100 mV–200 mV in Fig. 2. The outcome is explained by porous electrode theory and the distribution of Faradaic current throughout the cathode thickness. For the M–H and MHC mechanisms, the current distribution becomes more uniform throughout the cathode thickness, as j approaches jLim. On a local basis, the overpotential increases as needed to drive faster charge transfer rates. This also drives lithium ions to/from cathode regions closer to the current collector. Near the limiting current density, the transport resistance is lower than the kinetic resistance, and hence, increased transport flux dominates over any increase in the local current density. Ultimately, the Faradaic current distribution throughout the cathode will be the one that minimizes the total cathode overpotential. As the total current density approaches the limiting current value for a given mechanism, the increased current is therefore matched via small increases in the transport losses, rather than the large increases in the activation overpotentials predicted by Fig. 2.

This is demonstrated in Fig. 10, which shows the degree of lithiation non-uniformity for the BV and MHC P2D results. The degree of non-uniformity is here represented by the cathode particle surface lithium mole fraction for cathode particles at the electrolyte separator boundary (XLi,sep), relative to that at the current collector boundary (XLi,cc). The results shown in Fig. 10 are for a 40 μm cathode discharged at 30 mA cm−2, as in Fig. 8(b). At the beginning of discharge, all particles are uniformly delithiated. As cathode lithiation proceeds, we see that the particles at the separator are preferentially lithiated due to electrolyte transport limitations. This non-uniformity is present in both BV and MHC kinetics but is more pronounced for the BV predictions, where the charge transfer resistance is lower. As discharge proceeds and cathode particles near the separator approach full lithiation, the remaining lithiation current is more uniformly distributed. Compared to the discharge curve in Fig. 8(b), we see that this corresponds to an increasing difference in the MHC and BV potential profiles toward the very end of discharge.

FIG. 10.

Degree of lithiation non-uniformity in the cathode, calculated as the cathode particle surface lithium mole fraction at the electrolyte separator boundary (XLi,sep), relative to that at the current collector boundary (XLi,cc). Results show greater non-uniformity for Butler–Volmer kinetics relative to Marcus–Hush–Chidsey, reflecting higher charge transfer resistance in the M–H kinetics as the current density approaches the limiting value. Results are shown for a 40 μm cathode discharged at 30 mA cm−2, as in Fig. 8(b).

FIG. 10.

Degree of lithiation non-uniformity in the cathode, calculated as the cathode particle surface lithium mole fraction at the electrolyte separator boundary (XLi,sep), relative to that at the current collector boundary (XLi,cc). Results show greater non-uniformity for Butler–Volmer kinetics relative to Marcus–Hush–Chidsey, reflecting higher charge transfer resistance in the M–H kinetics as the current density approaches the limiting value. Results are shown for a 40 μm cathode discharged at 30 mA cm−2, as in Fig. 8(b).

Close modal

The results presented herein are, of course, for a single electrolyte composition (1M LiPF6 in ED:DEC) and cathode design (εelyte = 0.3). Changes to any of these properties (e.g., electrolyte salt concentration, cathode thickness, particle size, and microstructure) can decrease the transport resistance, relative to kinetic processes, leading to additional conditions where discrepancies between the MHC and BV kinetics are noteworthy. The existence of convenient, closed-form, computationally efficient MHC algorithms significantly reduces the barrier to the use of MHC kinetics in situations where they are most appropriate.

In this work, we used transient voltammetry data for lithium electrodeposition and stripping from the work of Boyle et al.10 to explore the compatibility of different kinetic theories to the data. The kinetic theories are compared in the context of developing better kinetic models for use in electrochemical performance simulations of lithium batteries. The widely used Butler–Volmer theory is found to be incompatible with the voltammetry data at high rates. We compare the use of the M–H kinetic model, proposed by Boyle et al., with the MHC model to explain the voltammetry data for LiPF6 in four solvents: PC, DEC, EC:DEC, and EC:DEC wt. 10% FEC. We find that both the M–H model and the MHC model fit the data with equal accuracy. The exchange current densities fit for each of the solvents are similar for the M–H and the MHC models.

We find that the reorganization energies predicted in this work for the M–H model match the reorganization energies predicted by Boyle et al. However, we find that the MHC model consistently predicts a (0.12 ± 0.1 eV) lower reorganization energy. We attribute this difference to the incorporation of the distribution of electrons at different energy levels within the MHC model, which the M–H model ignores. This is also consistent with other studies43 where it was found that the reorganization energy predicted by the MHC model was more accurate and that Marcus and M–H models required higher reorganization energies to fit Tafel curves. Given that the MHC-predicted reorganization energies are more accurate than Marcus or M–H predictions for other materials, we believe that this could hold for the lithium/solvent systems studied here as well. We also note that the M–H model proposed by Boyle et al. is a low overpotential approximation of the MHC model, and hence, the two models are indistinguishable at low overpotentials. At higher overpotentials, we observe significant deviations between the predictions of the M–H and MHC models. Thus, we propose the use of the MHC model implemented using a uniformly valid closed form approximation within battery models for describing the kinetics at lithium electrodes.

This recommendation is supported by P2D battery simulations of the extreme fast charge of a lithium metal anode paired with a LFP/C composite cathode. The observed differences in P2D predictions as a function of a kinetic mechanism are primarily due to kinetic limitations in the LFP cathode. While the kinetic mechanism disparities for Li are apparent only at current densities that exceed technologically relevant values, for LFP, the overpotentials predicted by the kinetic mechanisms diverge at C-rates of roughly 10 mA cm−2 to 30 mA cm−2. M–H kinetics predict that the battery cannot sustain charge–discharge at current densities greater than 1 mA cm−2 to 10 mA cm−2, in contrast to the available experimental evidence. This discrepancy is consistent with the theoretical underpinning of MH kinetics, which neglects the availability of valence-band electrons. The P2D results using BV and MHC kinetics, on the other hand, are quite similar to one another, except in kinetically limited regimes (i.e., at high current densities in thin cathodes), where MHC kinetics predict limiting current densities. Using an example eVTOL mission, we demonstrate that the MHC and BV models predict similar delivered power during the mission, except at the very beginning of takeoff and the very end of landing, where BV kinetics predict power densities roughly 1 mW cm−2 higher than MHC kinetics.

For galvanostatically controlled applications where the current densities do not exceed the MHC-predicted limiting current densities, BV and MHC kinetics are roughly equivalent. Much of this is due to a more uniform distribution of charge-transfer current through the cathode thickness, for MHC kinetics. However, the inability of the more commonly used BV and M–H kinetics to accurately predict limiting current densities cautions against their use in model-predicted design and control of batteries for high C-rate and high power applications. While the BV and M–H models offer the advantage of computational simplicity, closed-form approaches for MHC are available,11 and we demonstrate here that they can be readily implemented for validation against experimental data. To this end, we have published the model and thermo-kinetic modeling tools used here as open-source software to facilitate and accelerate the adoption of MHC kinetics more broadly throughout the field.46 

See the supplementary material for (i) additional visualizations of the fitting of kinetic models, (ii) a comparison between predictions for the reorganization energies for the M–H and MHC models, (iii) the use of the asymmetric Marcus–Hush model to fit the electrodeposition data, (iv) model parameters, and (v) model validation.

S.S. and D.K. contributed equally to this manuscript.

S.S. and V.V. gratefully acknowledge funding support from Aurora Flight Sciences. D.M.K. and S.C.D. gratefully acknowledge funding support from the National Science Foundation under Grant No. 1931584, Program Officer Raymond Adomaitis and from the Colorado Energy Research Collaboratory under Seed Grant No. 16-2018. S.C.D. and D.M.K. also thank Professor Wolfgang Bessler and Lutz Schiffer, Offenburg University of Applied Sciences, for insightful conversations about charge-transfer software implementations. The authors thank the Battery Modeling Webinar Series (BMWS) community for nucleating this collaboration and Battery Modeling Slack Workspace for help in managing this collaboration.

The data that support the findings of this study are openly available on Github at http://doi.org/github.com/coresresearch/BatCan/tree/MHC/li_ion.46 

1.
S.
Sripad
and
V.
Viswanathan
, “
Performance metrics required of next-generation batteries to make a practical electric semi truck
,”
ACS Energy Lett.
2
,
1669
1673
(
2017
).
2.
A.
Bills
,
S.
Sripad
,
W. L.
Fredericks
,
M.
Singh
, and
V.
Viswanathan
, “
Performance metrics required of next-generation batteries to electrify commercial aircraft
,”
ACS Energy Lett.
5
,
663
668
(
2020
).
3.
P.
Albertus
,
S.
Babinec
,
S.
Litzelman
, and
A.
Newman
, “
Status and challenges in enabling the lithium metal electrode for high-energy and low-cost rechargeable batteries
,”
Energy
3
,
16
21
(
2017
).
4.
A.
Mistry
and
V.
Srinivasan
, “
On our limited understanding of electrodeposition
,”
MRS Adv.
4
,
2843
2861
(
2019
).
5.
T.
Krauskopf
,
F. H.
Richter
,
W. G.
Zeier
, and
J.
Janek
, “
Physicochemical concepts of the lithium metal anode in solid-state batteries
,”
Chem. Rev.
120
,
7745
(
2020
).
6.
W. L.
Fredericks
,
S.
Sripad
,
G. C.
Bower
, and
V.
Viswanathan
, “
Performance metrics required of next-generation batteries to electrify vertical takeoff and landing (VTOL) aircraft
,”
ACS Energy Lett.
3
,
2989
2994
(
2018
).
7.
K. N.
Wood
,
E.
Kazyak
,
A. F.
Chadwick
,
K.-H.
Chen
,
J.-G.
Zhang
,
K.
Thornton
, and
N. P.
Dasgupta
, “
Dendrites and pits: Untangling the complex behavior of lithium metal anodes through operando video microscopy
,”
ACS Cent. Sci.
2
,
790
801
(
2016
).
8.
J.
Kasemchainan
,
S.
Zekoll
,
D.
Spencer Jolly
,
Z.
Ning
,
G. O.
Hartley
,
J.
Marrow
, and
P. G.
Bruce
, “
Critical stripping current leads to dendrite formation on plating in lithium anode solid electrolyte cells
,”
Nat. Mater.
18
,
1105
1111
(
2019
).
9.
M.
Doyle
,
T. F.
Fuller
, and
J.
Newman
, “
Modeling of galvanostatic charge and discharge of the lithium/polymer/insertion cell
,”
J. Electrochem. Soc.
140
,
1526
(
1993
).
10.
D. T.
Boyle
,
X.
Kong
,
A.
Pei
,
P. E.
Rudnicki
,
F.
Shi
,
W.
Huang
,
Z.
Bao
,
J.
Qin
, and
Y.
Cui
, “
Transient voltammetry with ultramicroelectrodes reveals the electron transfer kinetics of lithium metal anodes
,”
ACS Energy Lett.
5
,
701
709
(
2020
).
11.
Y.
Zeng
,
R. B.
Smith
,
P.
Bai
, and
M. Z.
Bazant
, “
Simple formula for Marcus–Hush–Chidsey kinetics
,”
J. Electroanal. Chem.
735
,
77
83
(
2014
).
12.
D. G.
Goodwin
,
R. L.
Speth
,
H. K.
Moffat
, and
B. W.
Weber
(
2018
). “
Cantera: An object-oriented software toolkit for chemical kinetics, thermodynamics, and transport processes
,” Zenodo, version 2.4.0, https://www.cantera.org.
13.
A. J.
Bard
,
L. R.
Faulkner
 et al,
Electrochemical Methods: Fundamentals and Applications
(
Wiley
,
New York
,
2001
).
14.
J. A. V.
Butler
, “
The mechanism of overvoltage and its relation to the combination of hydrogen atoms at metal electrodes
,”
Trans. Faraday Soc.
28
,
379
382
(
1932
).
15.
T.
Erdey-Gruz
and
M.
Volmer
, “
Zur theorie der wasserstoff Überspannung” (“The theory of hydrogen overvoltage
”),
Z. Phys. Chem.
150
,
203
213
(
1930
).
16.
R. A.
Marcus
, “
On the theory of oxidation-reduction reactions involving electron transfer. I
,”
J. Chem. Phys.
24
,
966
978
(
1956
).
17.
N. S.
Hush
, “
Adiabatic rate processes at electrodes. I. Energy-charge relationships
,”
J. Chem. Phys.
28
,
962
972
(
1958
).
18.
N. S.
Hush
, “
Electron transfer in retrospect and prospect: 1: Adiabatic electrode processes
,”
J. Electroanal. Chem.
460
,
5
29
(
1999
).
19.
C. E. D.
Chidsey
, “
Free energy and temperature dependence of electron transfer at the metal-electrolyte interface
,”
Science
251
,
919
922
(
1991
).
20.
J.
Newman
and
K. E.
Thomas-Alyea
,
Electrochemical Systems
(
John Wiley & Sons
,
2012
).
21.
V.
Ramadesigan
,
P. W. C.
Northrop
,
S.
De
,
S.
Santhanagopalan
,
R. D.
Braatz
, and
V. R.
Subramanian
, “
Modeling and simulation of lithium-ion batteries from a systems engineering perspective
,”
J. Electrochem. Soc.
159
,
R31
(
2012
).
22.
C.
Monroe
and
J.
Newman
, “
The impact of elastic deformation on deposition kinetics at lithium/polymer interfaces
,”
J. Electrochem. Soc.
152
,
A396
(
2005
).
23.
V.
Srinivasan
,
K.
Higa
,
P.
Barai
, and
Y.
Xie
, “
Computational modeling of morphology evolution in metal-based battery electrodes
,” in
Handbook of Materials Modeling: Methods: Theory and Modeling
, edited by
W.
Andreoni
and
S.
Yip
(
Springer International Publishing
,
Cham
,
2020
), pp.
1193
1219
.
24.
F.
Hao
,
A.
Verma
, and
P. P.
Mukherjee
, “
Mesoscale complexations in lithium electrodeposition
,”
ACS Appl. Mater. Interfaces
10
,
26320
26327
(
2018
).
25.
E.
Laborda
,
M. C.
Henstridge
,
C.
Batchelor-McAuley
, and
R. G.
Compton
, “
Asymmetric Marcus–Hush theory for voltammetry
,”
Chem. Soc. Rev.
42
,
4894
4905
(
2013
).
26.
J.
Savéant
and
D.
Tessier
, “
Convolution potential sweep voltammetry V. Determination of charge transfer kinetics deviating from the Butler-Volmer behaviour
,”
J. Electroanal. Chem.
65
,
57
66
(
1975
).
27.
R. A.
Marcus
and
N.
Sutin
, “
Electron transfers in chemistry and biology
,”
Biochim. Biophys. Acta
811
,
265
322
(
1985
).
28.
R. A.
Marcus
, “
On the theory of electron-transfer reactions. VI. Unified treatment for homogeneous and electrode reactions
,”
J. Chem. Phys.
43
,
679
701
(
1965
).
29.
M. C.
Henstridge
,
E.
Laborda
,
N. V.
Rees
, and
R. G.
Compton
, “
Marcus–Hush–Chidsey theory of electron transfer applied to voltammetry: A review
,”
Electrochim. Acta
84
,
12
20
(
2012
).
30.
R.
Nissim
,
C.
Batchelor-McAuley
,
M. C.
Henstridge
, and
R. G.
Compton
, “
Electrode kinetics at carbon electrodes and the density of electronic states
,”
Chem. Commun.
48
,
3294
3296
(
2012
).
31.
R.
Kurchin
and
V.
Viswanathan
, “
Marcus–Hush–Chidsey kinetics at electrode–electrolyte interfaces
,”
J. Chem. Phys.
153
,
134706
(
2020
).
32.
W. F.
Howard
and
R. M.
Spotnitz
, “
Theoretical evaluation of high-energy lithium metal phosphate cathode materials in Li-ion batteries
,”
J. Power Sources
165
,
887
891
(
2007
).
33.
D.
Wang
,
H.
Li
,
S.
Shi
,
X.
Huang
, and
L.
Chen
, “
Improving the rate performance of LiFePO4 by Fe-site doping
,”
Electrochim. Acta
50
,
2955
2958
(
2005
).
34.
Y.
Meng
,
Y.
Li
,
J.
Xia
,
Q.
Hu
,
X.
Ke
,
G.
Ren
, and
F.
Zhu
, “
F-doped LiFePO4@N/B/F-doped carbon as high performance cathode materials for Li-ion batteries
,”
Appl. Surf. Sci.
476
,
761
768
(
2019
).
35.
T. R.
Jow
,
M. B.
Marx
, and
J. L.
Allen
, “
Distinguishing Li+ charge transfer kinetics at NCA/electrolyte and graphite/electrolyte interfaces, and NCA/electrolyte and LFP/electrolyte interfaces in Li-ion cells
,”
J. Electrochem. Soc.
159
,
A604
A612
(
2012
).
36.
X.
Zhou
,
F.
Wang
,
Y.
Zhu
, and
Z.
Liu
, “
Graphene modified LiFePO4 cathode materials for high power lithium ion batteries
,”
J. Mater. Chem.
21
,
3353
3358
(
2011
).
37.
H.
Li
,
L.
Peng
,
D.
Wu
,
J.
Wu
,
Y.-J.
Zhu
, and
X.
Hu
, “
Ultrahigh-capacity and fire-resistant LiFePO4-based composite cathodes for advanced lithium-ion batteries
,”
Adv. Energy Mater.
9
,
1802930
(
2019
).
38.
A. M.
Colclasure
and
R. J.
Kee
, “
Thermodynamically consistent modeling of elementary electrochemistry in lithium-ion batteries
,”
Electrochim. Acta
55
,
8960
8973
(
2010
).
39.
J.
Chiew
,
C. S.
Chin
,
W. D.
Toh
,
Z.
Gao
,
J.
Jia
, and
C. Z.
Zhang
, “
A pseudo three-dimensional electrochemical thermal model of a cylindrical LiFePO4/graphite battery
,”
Appl. Therm. Eng.
147
,
450
463
(
2019
).
40.
M. D.
Murbach
and
D. T.
Schwartz
, “
Extending newman’s pseudo-two-dimensional lithium-ion battery impedance simulation approach to include the nonlinear harmonic response
,”
J. Electrochem. Soc.
164
,
E3311
E3320
(
2017
).
41.
C.
Kupper
and
W. G.
Bessler
, “
Multi-scale thermo-electrochemical modeling of performance and aging of a LiFePO4/graphite lithium-ion cell
,”
J. Electrochem. Soc.
164
,
A304
A320
(
2017
).
42.
P.
Lu
,
C.
Li
,
E. W.
Schneider
, and
S. J.
Harris
, “
Chemistry, impedance, and morphology evolution in solid electrolyte interphase films during formation in lithium ion batteries
,”
J. Phys. Chem. C
118
,
896
903
(
2014
).
43.
P.
Bai
and
M. Z.
Bazant
, “
Charge transfer kinetics at the solid–solid interface in porous electrodes
,”
Nat. Commun.
5
,
3585
(
2014
).
44.
Y.-C.
Chang
,
C.-T.
Peng
, and
I.-M.
Hung
, “
Effects of particle size and carbon coating on electrochemical properties of LiFePO4/C prepared by hydrothermal method
,”
J. Mater. Sci.
49
,
6907
6916
(
2014
).
45.
J.
Tang
,
X.
Zhong
,
H.
Li
,
Y.
Li
,
F.
Pan
, and
B.
Xu
, “
In-situ and selectively laser reduced graphene oxide sheets as excellent conductive additive for high rate capability LiFePO4 lithium ion batteries
,”
J. Power Sources
412
,
677
682
(
2019
).
46.
See https://github.com/coresresearch/BatCan/tree/MHC/li_ion for Batcan: Battery modeling with Cantera.
47.
48.
D. A.
Cogswell
, “
Quantitative phase-field modeling of dendritic electrodeposition
,”
Phys. Rev. E
92
,
011301
(
2015
).
49.
W. R.
Fawcett
, “
Potential dependence of the elementary steps in the kinetics of electrode reactions involving amalgam formation
,”
J. Phys. Chem.
93
,
2675
2682
(
1989
).
50.
Y.
Zeng
,
P.
Bai
,
R. B.
Smith
, and
M. Z.
Bazant
, “
Simple formula for asymmetric Marcus–Hush kinetics
,”
J. Electroanal. Chem.
748
,
52
57
(
2015
).

Supplementary Material