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 LiFePO_{4} 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.

## I. INTRODUCTION

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.

## II. METHODS

### A. Electron transfer kinetics

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 Marcus^{16} in the 1950s, kinetics of electrode processes examined by Hush,^{17,18} and, more recently, heterogeneous electron transfer explained using the Marcus–Hush–Chidsey^{19} 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:

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

where k_{red} and k_{ox} 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 Tessier^{26} 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

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

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 Tessier^{26} 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,

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 = j_{0}·(k_{red} − k_{ox}), where j_{0} is the exchange current density. The net current for the BV model can be obtained using Eqs. (2) and (3),

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 j^{M–H} would have its exchange current density $j0M\u2212H$ 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,

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 be^{11}

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 LiPF_{6} 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.

### B. Pseudo-2D fast-charge model

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 LiFePO_{4}–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.

Despite the lower gravimetric capacity and electrical conductivity of LiFePO_{4}, it has several benefits, such as lower material cost (6.3 Wh US$^{−1}) than other standard cathode materials^{32} 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 interface^{35} 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 LiFePO_{4} phase,

X

_{Li,ca}is the intercalation fraction in the cathode (−) andΦ

_{ca}is the electric potential of the cathode (V).

In the electrolyte phase,

X

_{k,elyte}is the mole fraction of species*k*in the electrolyte (kmol_{k}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 j_{ext}. 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 D_{Li,ca}, the governing equation for the intercalated lithium mole fraction is derived from conservation of mass and elements,

where C_{ca} is the total molar concentration of the cathode phase (mol m^{−3}), and the term $s\u0307Li$ is a molar source term due to the charge transfer reaction at the LFP particle surfaces (mol m^{−2} s^{−1}). A_{int} is the specific surface area of the electrode/electrolyte interface ($mint2\u2009m\u22123$). $\epsilon $_{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,

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, $s\u0307LiAint$ 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,

where C_{elyte} is the electrolyte molar concentration (mol m^{−3}, assumed to be constant for the incompressible electrolyte), $\epsilon $_{elyte} is the electrolyte volume fraction, *N*_{k, elyte} is the molar flux (mol_{k} m^{−2} s^{−1}), calculated according to dilute solution theory,^{41} and $s\u0307k,\u2009elyte$ is the net production rate of electrolyte species k ($molk\u2009mint\u22122\u2009s\u22121$) due to heterogeneous reactions at the electrode/electrolyte interface. Within the electrolyte separator, the term $s\u0307k,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

and

In these equations, j_{dl} is the double-layer current per unit volume at the electrode/electrolyte interface (A m^{−3}), C_{dl,ed} is the double-layer capacitance (F m^{−2}), and j_{io} and j_{el} 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),

where j_{Far} is the Faradaic current at the electrode/electrolyte interface due to charge transfer reactions, per unit volume (A m^{−3}). The interfacial currents j_{Far} and j_{dl} are both defined such that positive j represents positive charge entering the electrode and leaving the electrolyte. Therefore, positive j_{dl} 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 models^{38,39} assume that the double layer responds infinitely quickly such that j_{dl} 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 j_{Far} = −∇j_{el}. 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Ω m^{2} is used here^{42}). 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 j_{Far} = j_{SEI}. 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–LiFePO_{4}/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 literature^{45} 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}

## III. RESULTS

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, j_{0}’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.

### A. Electron transfer kinetics

The best fit exchange current densities for the M–H model ($j0M\u2212H$) 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 $j0M\u2212H$ 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.

Solvent . | $ j o M \u2212 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 \u2212 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.

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 Bazant^{43} 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 Bazant^{43} 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.

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 studies^{48,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/cm^{2} 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/cm^{2}. 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/cm^{2}. Similarly, for lower exchange current density solvents such as PC and DEC, the deviation is about 50 mA/cm^{2}. 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 approximation^{11} 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.

### B. Pseudo-2D fast-charge model

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.

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 j_{ext} and the total j_{Far} decreases, and the slope of the voltage profile also decreases. However, once the overpotentials reach the point corresponding to j_{Far, max} (roughly *η* = 0.35 V, in Fig. 2), further charging of the double layer decreases j_{Far}, 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 j_{Far} goes to zero with further increases in ΔΦ_{dl}, the double layer current is essentially constant at j_{dl} = j_{ext}, 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).

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 ≥ j_{Lim} 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 $gLFP\u22121$ at 20 to 30 mA cm^{−2}.

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 < j_{Lim}, 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.

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 j_{Lim} (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 j_{Lim}. 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 (*X*_{Li,sep}), relative to that at the current collector boundary (*X*_{Li,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.

The results presented herein are, of course, for a single electrolyte composition (1M LiPF_{6} in ED:DEC) and cathode design ($\epsilon $_{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.

## IV. CONCLUSIONS

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 LiPF_{6} 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 studies^{43} 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}

## SUPPLEMENTARY MATERIAL

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.

## AUTHORS’ CONTRIBUTIONS

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

## ACKNOWLEDGMENTS

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.

## DATA AVAILABILITY

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}

## REFERENCES

_{4}by Fe-site doping

_{4}@N/B/F-doped carbon as high performance cathode materials for Li-ion batteries

^{+}charge transfer kinetics at NCA/electrolyte and graphite/electrolyte interfaces, and NCA/electrolyte and LFP/electrolyte interfaces in Li-ion cells

_{4}cathode materials for high power lithium ion batteries

_{4}-based composite cathodes for advanced lithium-ion batteries

_{4}/graphite battery

_{4}/graphite lithium-ion cell

_{4}/C prepared by hydrothermal method

_{4}lithium ion batteries