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.
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 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:
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 kred and kox are the reaction rates for reduction and oxidation, respectively, 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 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
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 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,
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),
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 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 be11
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 and 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 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.
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,
where Cca is the total molar concentration of the cathode phase (mol m−3), and the term 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 (). 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, 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 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 is the net production rate of electrolyte species k () due to heterogeneous reactions at the electrode/electrolyte interface. Within the electrolyte separator, the term 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, 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),
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
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, 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.
A. Electron transfer kinetics
The best fit exchange current densities for the M–H model () and the MHC model (), along with the reported linear fit values () 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 values are seen to be marginally higher than and . 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.
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 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.
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.
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 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).
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 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 < 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.
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.
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.
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 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
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