Statistical Physics-Based Model of Solid Electrolyte Interphase Growth in Lithium Ion Batteries

The article presents a statistical physics-based model for the growth of the solid electrolyte interphase (SEI) in the negative electrode of lithium ion batteries. During battery operation, the SEI thickness grows by the reaction between lithium ions, electrons and solvent species on the surface of active particles at the negative electrode. The growth of the SEI layer causes a loss of lithium ions that induces capacity fade. In addition, it increases the ion transport resistance and decreases the total porosity. Our model employs a population balance formalism based on the Fokker-Planck Equation to describe the propagation of the particle density distribution function in the electrode. Structure-transforming processes at the level of individual particles are accounted for by using a set of kinetic and transport equations. These processes alter the particle density distribution function, and cause changes in battery performance. A parametric study of the model singles out the first moment of the initial SEI thickness distribution as the determining factor in predicting the capacity fade. The model-based treatment of experimental data allows analyzing processes that control SEI growth and extracting their controlling parameters. © The Author(s) 2017. Published by ECS. This is an open access article distributed under the terms of the Creative Commons Attribution 4.0 License (CC BY, http://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse of the work in any medium, provided the original work is properly cited. [DOI: 10.1149/2.1581706jes] All rights reserved.

Lithium ion batteries (LIBs) are highly touted energy storage devices for portable electronics and electric vehicles. 1 The LIB system consists of a negative electrode, a positive electrode, a separator, an electrolyte, and two current collectors. The most commonly used electrolytes are comprised of lithium salts, such as LiPF 6 in a solution of ethylene carbonate (EC) and dimethyl carbonate (DMC). 1 The electrodes consist of randomly distributed and interconnected particles of active material, which store and release the lithium ions.
Aging and degradation of LIBs have become major concerns for the operation of electric vehicles (EVs), which must fulfill exacting requirements in terms of durability, cyclability and overall lifetime. 2 Battery aging is linked to the (electro-)chemical and mechanical degradation of the electrodes and the electrolyte. 2 Three main degradation mechanisms prevail at the particle level during the cycling of LIBs: (1) growth of the solid electrolyte interphase (SEI) at the negative electrode, 2,3 (2) formation of cracks in the SEI layer at the negative electrode, 4 and (3) dissolution and isolation of nanocrystalline particles at the positive electrode. 5 These processes lead to capacity fade and power loss.
At the beginning of the battery life, particularly during the first cycle, the electrolyte undergoes reduction at the electrode/electrolyte interface, because the negative electrode operates at potentials that are outside of the electrochemical stability window of electrolyte components. This reduction is accompanied by the irreversible consumption of lithium ions. 6 This process forms a passivating surface layer known as the solid electrolyte interphase. The SEI layer grows further during charging, thereby reducing the amount of active lithium ions. Moreover, this layer penetrates into the electrode pores, reducing the overall porosity and decreasing ion access to the active surface area of the electrode. 2 Understanding the mechanism of SEI thickness growth and the resulting composition of the layer is a vital topic of LIB research, because of its great practical significance and the intricate interplay of underlying processes. 2,[6][7][8][9][10][11][12][13][14] The growth of the SEI involves a complex reaction network 15,16 that is sensitive to operating conditions as well as the battery cycling protocol. A mechanistic understanding of the processes involved in SEI formation and growth is vital for designing LIB systems with high performance and long cycle life. The topic therefore garners significant interest in the experimental and theoretical research community.
Recent efforts in modeling of the electrochemical processes during cycling of LIBs include particle level models 17,18 and porous electrode theory. 16 In 1976, Bennion developed the first model of the SEI for lithium metal electrodes. 19 Peled further derived a parabolic growth law for the SEI considering the electron transfer as the rate-determining step in the growth of the SEI layer. 20 Ploehn et al. developed a model for SEI growth, in which the diffusion of solvent through the SEI is the rate-determining step. 21 Interestingly, this assumption also leads to a parabolic growth law. Newman et al. developed a model to predict the SEI growth based on porous electrode theory, including the detailed chemistry of SEI formation. 16 This model was the first attempt to link chemical reactions leading to SEI formation to irreversible capacity loss. However, it comprises a large number of model parameters. Xie et al. expanded Newman's model to incorporate thermal effects. 22 They found the battery skin temperature to significantly affect the growth of the SEI. Several subsequent works are based on Newman's model. Their common feature is that they couple a one dimensional model of charge and mass transport through the electrode to a one dimensional model of spherical diffusion inside active particles, resulting in a pseudo two-dimensional porous electrode (P2D) model. [23][24][25][26] Under low to moderate working conditions, the P2D model can be reduced to a single particle (SP) model. 27,28 In this approach, the electrode is represented by identical spherical particles, whose accumulated surface area is equivalent to the total active area of the solid phase in the porous electrode. This model assumes that the mass and charge transport resistances in solution can be neglected. Moreover, the faradaic current across the porous electrode exhibits a linear profile, which is valid in the limit of small applied currents in thin and highly conductive electrodes. 18 The modeling work presented in this article strives to establish relations among structural properties and electrochemical performance of the battery in order to rationalize the influence of SEI growth upon them. The central component of the statistical physics-based modeling framework is the Fokker-Planck Equation (FPE). The porous battery electrode is treated as a statistical distribution of interconnected particles. The FPE governs the temporal evolution of this distribution. The developed formalism allows leading causes of structural degradation, driven by (electro-)chemical and mechanical stressors, to be incorporated. Capacity and power fade as well as the battery cycle life can be analyzed in dependence of the initial structure, the external conditions and the cycling protocol applied, to predict the propagation of the particle density distribution function due to SEI growth.
Eikerling and coworkers introduced a statistical modeling framework to study the degradation of platinum nanoparticles in the cathode catalyst layer of a polymer electrolyte fuel cell. [29][30][31][32] In the area of LIB research, Dreyer et al. employed a population balance approach to study the phase-change behavior of randomly dispersed interconnected particles during lithium storage. 33 Röder et al. used a similar approach to study the influence of the particle size distribution on LIB performance. 34 The present work aims to investigate the evolution of the SEI thickness distribution during cycling of the battery under different aging conditions. The following section introduces the statistical electrode model and it describes how the SP model for SEI growth is incorporated. Details of the numerical solution and the fitting strategy for the treatment of experimental data are given; a discussion of possible discretization schemes is provided in the Appendix. The section Results and Discussion rationalizes the dependence of capacity fade on the average SEI thickness as well as the width and tail of the SEI thickness distribution. The capabilities of the model are demonstrated by comparison to experimental data for capacity fade at varying temperature and charging rate.

Model Development
In our model, the electrodes are assumed to consist of random distributions of active particles. Particles at both electrodes are connected by inert binders. Further assumptions are that (1) the concentrations of solvent and Li ions in the solution phase are uniform, (2) the solution phase potential is uniform, (3) the electric potential in the solid phase of active particles is uniform; and (4) the faradaic current is distributed uniformly on the particle surfaces. Given these assumptions, particles charge and discharge uniformly throughout the electrode at all C-rates. The possibility of non-uniform current distribution or a particle-by-particle charging mechanism as discussed by Li et al. 35 will be considered in a forthcoming work.
The transport of solvent molecules inside the SEI layer proceeds via diffusion, as shown in Fig. 1. The material balance for the solvent concentration inside of the SEI layer is The boundary conditions for this equation are and where D s is the diffusion coefficient for solvent molecules in the SEI layer and C s,b is the solvent bulk concentration. The diffusion of lithium ions in an active solid particle is governed by Fick's second law in spherical coordinates, where C denotes the concentration of lithium ions in the particle, D i is the diffusion coefficient of Li ions in the solid and r is the radial coordinate. The boundary conditions are For the sake of simplicity, we assume that the SEI layer grows in a reaction between Li ions, electrons and solvent molecules, where Solv represents a solvent molecule such as ethylene carbonate, dimethyl carbonate, or diethyl carbonate and P stands for the reaction product. The evolution of the SEI thickness is described by the Faraday law, The SEI growth reaction is assumed to be irreversible and to follow a Tafel kinetics, In this equation, k SE I is the reaction rate constant for SEI growth, C s is the solvent concentration at the particle surface and α is the electronic charge transfer coefficient. The overpotential of SEI growth is In Eq. 10, ϕ M and ϕ E are the potentials in metal phase and solution phase, respectively, and κ SE I is the ionic conductivity of the SEI. U re f SE I is the equilibrium potential for the SEI growth, for which a value of 0.4 V vs. Li electrode has been reported. 17 J t in Eq. 10 is the total current density, The electrochemical reaction for intercalation/deintercalation of Li ions at the solid/solution interface can be written as where θ s represents a vacant site on the particle surface. 18 The rate of this reaction is given by the Butler-Volmer Equation 16, where k int is the intercalation reaction rate constant, C max the maximum solid phase concentration, C sur f the Li ion concentration at the surface of active particles and C e the concentration of Li ions in the electrolyte, which is assumed to be constant. The overpotential for the intercalation reaction at the negative electrode, η N int , is given by Here, U N is the open circuit voltage of the negative electrode, which was extracted from experimental data of Liu et al. 36 by Pinson et al. 37

Population Balance Model
In this work, we employ a population balance model (PBM) based on the Fokker-Planck Equation. This approach provides a predictive framework to explore the structural evolution of a statistical distribution of solid particles, droplets, or bubbles. 38 Within this framework, a random distribution of objects propagates according to the kinetics of processes that transform the properties of individual objects. The rate of these transformations of objects depends on local conditions inside the medium. 39 Examples for systems described successfully by this formalism include catalyst layers of PEM fuel cells, in which the catalyst particle radius distribution changes as a result of degradation, [29][30][31][32]40 biological cells in biomedical systems, e.g. in cancer research, 41 or monomer droplets during suspension polymerization of vinyl chloride. 42 In the mathematical formalism, the electrode is represented by a particle density distribution f ( X , t), which is the particle number per unit volume of the electrode in the parameter space spanned by the vector variable X . The temporal evolution of f ( X , t) is governed by the FPE, 43 where X represents a set of quantities that describe the state of individual particles in the distribution, including for instance particle radius, SEI thickness or charging state. The drift term d X dt accounts for the continuous transformation of these particle properties. Source/sink terms, included in k Q k , represent processes that delete or add particles from or to the distribution.
The present work specifically focuses on the growth of the SEI layer thickness as the structure transforming process at the single particle level. We therefore replace X by the SEI layer thickness, δ. Moreover, it is assumed that the total number of particles is constant, which implies that source/sink terms are neglected. With these assumptions, we obtain a simplified FPE, where dδ dt represents the SEI thickness change. Given an initial function f (δ, 0), Eq. 16 allows the propagation of the SEI thickness and of dependent properties to be computed. The accuracy of the calculation depends on the knowledge of f (δ, 0) and on the inclusion and suitable parameterization of particle-level processes. We consider f (δ, 0) as the distribution of the SEI thickness formed during the production of the battery; this distribution depends on the fabrication approach and conditions as well as the conditioning parameters of the battery. For the present work, we determined the initial distribution using data from SEM image analysis of a LIB negative electrode, reported in Ref. 11. These data were interpolated by an exponentially modified Gaussian (EMG) function, given by [17] with h = 0.9, μ = 32.0, σ = 8.0 and ρ = 7.5 obtained from fitting.
The solution of Eq. 16 allows for the calculation of transiently changing moments of the distribution, viz. average SEI thickness,δ, and capacity fade, Q loss , as given bȳ [20] where M SE I and ρ SE I are the molar mass and the density of the SEI film and Q 0 is the initial capacity of the battery.

Model Solution and Data Fitting
The system of equations is solved in MATLAB, using a central difference discretization scheme to transform the PDEs into a system of ordinary differential equations (ODEs). All ODEs are solved simultaneously employing the ode15s solver in MATLAB, which is a variable-step, variable-order solver based on the numerical differentiation formulae of order 5. Discretization of the PDEs is explained in the Appendix. Instead of time, we consider the number of cycles, N , as the independent transient variable, considering kinetic model parameters as effective parameters, i.e. average values over a single cycle.
The model was compared to capacity fade data from Wang et al. for varying T and charging rate. We adopted the same method for the capacity fade calculation as Wang et al., 44 i.e., the capacity measured during the first cycle is used as Q 0 . For the capacity fade simulations at different temperatures, a fitting procedure was implemented in MAT-LAB to find k SE I and D Solv for one set of operating conditions (15 • C and C/2). A charging rate of 1C corresponds to a current of 2 A. 44 At the other two temperatures (45 • C and 60 • C), only k SE I was fitted. For the simulation of capacity fade at different charging rates, the fitting procedure was used to obtain k SE I and D s at 20 • C and 2C. At the other two charging rates (6C, 10C), only k SE I was fitted. The list of parameters used in the simulations is presented in Table I. Parameters are based on the experimental work of Wang et al. 44 and they were extracted by Prada et al. 48

Results and Discussion
Evolution of SEI layer.-We first evaluated the impact of the SEI thickness growth on the single particle properties. From the initial SEI distribution (Eq. 17), we chose three sample particles with initial SEI thicknesses of 56 nm, 83 nm and 116 nm to simulate their SEI growth at 15 • C and under a C/2 charging rate, as shown in Fig. 2. Initially, the SEI layer grows quickly. At higher cycle numbers, the rate of growth slows down due to the increasing ionic and diffusion resistance of the SEI layer. Since SEI thickness growth is slower for particles with initially thicker SEI layer, the three curves converge.
The evolution of the SEI thickness distribution during LIB cycling is shown in Fig. 3. The overarching trend is a shift of the particle density distribution to larger average SEI thickness. Since the SEI thickness grows faster on particles with smaller initial SEI thickness, the SEI distribution narrows with N , with a corresponding increase of the height of the peak.
Effect of initial SEI distribution.-As explained in the modeling section, we use an EMG function (Eq. 17) to describe the initial SEI distribution. Fig. 4 shows the EMG function along with the SEM image analysis data. Three features determine the SEI distribution: (1) the average SEI thickness,δ, (2) the width, (3) and the tail of the distribution. These features can be influenced during the manufacturing process and initial conditioning of the battery. We have studied the effect of each factor on the shape of the SEI distribution and on the capacity fade.
For comparative analyses of the impact of these characteristics, we created four initial SEI distributions, as shown in Fig. 5a: an EMG distribution as a reference case withδ = 40 nm; two EMG distributions with largerδ, which were created by increasing μ from 32 to 42 and 52; an EMG distribution with the sameδ as the reference case but different values of σ and ρ. The capacity fade incurred by the SEI formation that occurs during the first charge-discharge cycle of the freshly made battery and during SEI growth over 1000 subsequent cycles was simulated at 15 • C and under a C/2 charging rate. The results are shown in Fig. 5b.   Figure 2. Growth of the SEI layer for a single particle at 15 • C and C/2 charging rate. We found thatδ exerts the most significant impact on Q loss , whereas the width and the tail of the particle distribution did not exert any discernible effect. This can be seen in Fig. 5b. The curves of Q loss for the two EMG distributions withδ = 40 nm but different σ and ρ lie on top of each other. Fig. 6 shows the result of fitting the model to experimental data for Q loss of Wang et al. 44 at T = 15 • C, 45 • C, and 60 • C in Fig. 6a and at charging rates of 2C, 6C and 10C in Fig. 6b. The fitted curves were obtained by variation of k SE I and D s . Only the capacity fade caused by SEI growth was included for the comparison with experimental data. Higher T accelerates SEI growth due to the enhanced electrochemical reaction rate. The obtained model fits are in very good agreement with experimental data. From the Arrhenius plot in Fig. 7a and Fig. 7b we calculated an activation energy of 0.45 eV for the SEI growth reaction and 0.02 eV for solvent diffusion through the SEI layer. The activation energy of the SEI growth reaction in Eq. 7 is highly dependent on the surface properties of the negative electrode material and more importantly on the solvent properties. 45 Different values have been reported in the literature, ranging from 0.34 eV to 1.5 eV. [45][46][47] The scatter of activation energy values can be attributed to variations in cycling conditions as well as differences in materials structure and composition. For instance, we do not know what type of solvent was used by Wang et al. 44     Fig. 8a shows the change in the SEI thickness distribution after 100 cycles for the three curves in Fig. 6a. When T increases, the distribution shifts to higherδ and the distributions becomes narrower. Fig. 6b reveals excellent agreement of the model fits with experimental data of Wang et al. 44 for Q loss at the different C-rates. The obtained values for k SE I are 1.18 × 10 −17 m/s at 2C, 1.47 × 10 −17 m/s at 6C, and 3.59 × 10 −17 m/s at 10C. The effect of the charging rate on the change in the SEI thickness distribution after 100 cycles at 15 • C is shown in Fig. 8b. By increasing the C-rate, the SEI distribution is shifted to largerδ and the narrowing of the distribution, described above, becomes more pronounced.

Conclusions
This article presented a statistical physics-based model to study the impact of SEI growth on the structure and performance of the negative electrode of lithium ion batteries. The two-scale approach employs a formalism based on the Fokker-Planck Equation to relate the SEI thickness growth at the single particle level to the changes in structural properties at the electrode level. The model allows exploring the correlations of structural changes that occur during charge-discharge cycling and electrode performance.
A parametric analysis of the model reveals that the first moment of the initial SEI distribution, i.e., the average SEI thickness, is the determining factor in predicting the capacity fade. The shape and width of the initial SEI distribution do not exert any discernible impact on capacity fade. During battery cycling, the SEI thickness distribution of particle-based electrodes converges toward a narrowly shaped distribution with uniform SEI thickness.
Comparison of the model to experimental data showed convincing agreement and it allows for the rate of SEI growth and the solvent diffusion coefficient to be obtained. Further analysis using Arrhenius plots provides the activation energies of these processes. Given a set of initial structural and kinetic parameters for a specific battery electrode, the model will be able to predict the evolution of electrode structure and performance upon exposure of the battery to varying life cycle protocols. The versatile modeling framework will allow incorporating further structure-transforming processes in order to systematically improve the accuracy of that prediction.