Fractional Differential Generalization of the Single Particle Model of a Lithium-Ion Cell

The effect of anomalous diffusion of lithium on the discharge curves and impedance spectra of lithium-ion batteries (LIB) is studied within the fractional differential generalization of the single-particle model. The distribution of lithium ions in electrolyte and electrode particles is expressed through the Mittag–Leffler function and the Lévy stable density. Using the new model, we generalize the equivalent circuit of LIB. The slope of the low-frequency rectilinear part of the Nyquist diagram does not always unambiguously determine the subdiffusion index and can be either larger or smaller than the slope corresponding to normal diffusion. The new aspect of capacity degradation related to a change in the type of ion diffusion in LIB components is discussed.


Introduction
Mathematical modeling of lithium-ion batteries (LIB) is necessary both to optimize the design and operation modes and to test the state of working devices: to assess state-of-health, state-of-charge, and predict performance [1,2]. The charging and discharging of LIBs are followed by complex processes including reaction at the electrode/electrolyte interface, lithium transport in solid electrode particles, the formation of an electrical double layer, the growth of a solid-electrolyte interface (SEI), phase transformations, and redistribution of stresses during lithiation and delithiation, among others.
Recently, LIB models based on fractional-order impedances have been proposed [3][4][5]. These models are justified by the anomalous diffusion of ions in a disordered medium of a percolative type. In [6], the authors used the phenomenological model accounting for anomalous ion diffusion in a LiCoO 2 cathode to describe the low-frequency part of the Nyquist diagram that has a steeper slope than predicted by normal diffusion. Using the subdiffusion generalization [7] of the electrochemical impedance, the authors of [6] showed that the anomalous-diffusion model describes the impedance of LiCoO 2 |C cells well at normal operating potentials and temperatures. They noted that the contribution of anomalous diffusion decreases for higher temperatures. Some experimental studies (see, e.g., [8][9][10]) using nuclear magnetic resonance (NMR) clearly showed that the diffusion of ions and molecular complexes in some media is characterized by features typical of anomalous diffusion. The mathematical basis of anomalous self-similar diffusion is formed by equations with fractional-order derivatives [11][12][13]. Sabatier et al. [14,15] argued that fractional differentiation helps to describe some components of a model with a smaller number of parameters related to electrochemical quantities. Such fractional models with a reduced number of fitting parameters can be used for the development of state-of-charge (SOC) and state-of-health (SOH) estimators. In [14,15], fractional models were used as simplifications of an integer-order electrochemical model, and the authors did not consider the anomalous diffusive behavior of ion transport in electrolyte and electrode particles.
Some solid electrolytes used in LIB, for example, compound Li 7 La 3 Zr 2 O 12 (LLZO), have a polycrystalline structure with grain boundaries (GBs). As is known, GBs can accelerate or suppress the diffusion of atoms and ions in polycrystalline materials (see [16]). Many approaches to the description of grain boundary diffusion are based on the classical model proposed by Fisher in 1951. They use a system of normal-diffusion equations in a heterogeneous medium. A direct relationship between the fractional differential model of anomalous diffusion and the Fisher model and its modifications without additional assumptions was derived in [17] (see also [18]).
Many models of ion transport in a porous electrode are based on the effective medium approximation, where two parameters are introduced: the volume fraction of space available for transfer and the tortuosity of flow paths. With the use of fractal geometry and fractional calculus, Huang [19] proposed a new model to describe the mass transfer and reactions in LIB electrodes having a fragmented fractal structure.
Despite the observed evidence of anomalous diffusion and the success of fractional circuit models, there are no microscopic models that jointly describe anomalous diffusion and percolation of lithium ions in porous electrodes, SEI, and electrode particles. In this work, we present an attempt to make such a generalization using the fractional modification of the single-particle (SP) model.

Fractional Differential Diffusion Equation
Anomalous diffusion is often characterized by the power law expansion of the diffusion packet, ∆(t) ∝ t α/2 . The classical diffusion (α = 1) implies Brownian motion of particles and the Fickian relation between current and concentration. When 0 < α < 1, the diffusion packet width increases with time more slowly than in the normal case (subdiffusion). If α > 1, it increases more quickly (superdiffusion). The mathematical basis of anomalous self-similar diffusion is formed by equations with fractional-order derivatives [11][12][13]. It should be noted, however, that linear growth of the diffusion packet width does not always correspond to ordinary diffusion (see, e.g., [20]) and can be observed for processes with "diffusing diffusivities" characterized by non-Gaussian distributions.
The dispersive diffusion (when delay is associated with localization times) can be described by the following integro-differential equation (see the reviews [11,13,21]): where Q(t) is a memory kernel and c(r, t) is a particle concentration. In the subdiffusive regime, the requirement of self-similarity of the Green's functions leads to the fractional-order diffusion equation [11]: Here: is the fractional Riemann-Liouville derivative of order ν [22], Γ(ν) is the gamma function, and K is a subdiffusion coefficient. Fundamental solutions to this equation can be found in [11]. Fractional differential Equation (2) can be obtained as the asymptotic equation describing the continuous time random walk (CTRW) [13] with "heavy-tailed" distributions of random localization times τ: The power law distribution of localization times can be obtained, for example, in the random activation energy model. Suppose that ion hopping frequency is defined by a semiclassical expression W = Ae −∆E/kT , and the conditional distribution of waiting time corresponding to a given activation energy ∆E = ε is exponential: P{T > t|ε} = e −W(ε)t . If the activation energy is a random variable with the following distribution density p(ε) = ε −1 0 e −ε/ε 0 , then averaging over activation energies leads to the waiting time distribution of the power law type: with α = kT/ε 0 . Here, ε 0 is a parameter characterizing the dispersion of activation energy. Under condition kT < ε 0 , ion transport is subdiffusive.
Equation (2) can be justified within other transport models (see, e.g., [23]). In electrodes based on porous and nanostructured materials, there are many alternative transport paths, and ion transport is often percolative. To study diffusion on percolation clusters, the comb model [24] is often used. In the simplest version [24,25], the order of the fractional derivative in Equation (2) is equal to 1/2. In more complex versions, for example, in the hierarchical comb model [26] or in the model with random finger lengths [27], other values of the subdiffusion parameter can be obtained.

Subdiffusion Model of the Lithium-Ion Cell
To assess the effect of anomalous diffusion on the discharge characteristics of LIB and impedance spectra, we consider the subdiffusive generalization of the single-particle (SP) model, which is a simplification of the pseudo-two-dimensional (P2D) model ( Figure 1). The original P2D-model was developed in [1,28,29]. The segment [0, L] on the x axis is considered to be divided into three sections corresponding to two electrodes and a separator. Newman and Tobias [30] showed that the current distribution in the porous electrode due to the electrochemical reaction is fairly uniform.
The assumption about homogeneous current across the electrode thickness is a key point of the model. In the standard SP model, electrodes are modeled as two spherical particles, and the variations in electrolyte concentration and in potential are usually ignored. Although, there are generalizations (see, e.g., [31,32]), where these variations are taken into account. Here, we also consider changes in concentration and potential along the x -axis. Spread in particle sizes can be taken into account within the so-called multiple particle model (generalization of the SP model). Meyers et al. [33], combining the impedance of a spherical electrode particle with an electrolyte impedance, obtained an analytical expression for the impedance of a porous electrode under the assumption about the ohmic character of electrolyte impedance. In [34], the authors proposed the impedance model describing diffusion within the solid electrode particles, charge transfer, charging of the double layer, formation of a solid-electrolyte interface, and the concentration gradients inside the electrolyte. Guo et al. [31] studied the energy balance within the SP model and described the thermal behavior of LIB in the process of galvanostatic charge-discharge modes. Note that existing versions of the P2D and SP models do not take the anomalous diffusive propagation of ions in porous electrodes, the solid-electrolyte interface, and amorphous phases into account. Variations of LIB discharge curves calculated within the SP model under normal diffusion conditions are demonstrated in Figure 2. The parameters and formula of the model are given below.     Table 1, subdiffusion indexes are taken equal to one.
We consider the general case of subdiffusive motion of lithium in spherical electrode particles and electrolyte. The fractional-order equation of isotropic subdiffusion in spherical coordinates has the form: where j indicates one of the electrodes (j = p for the cathode and j = n for the anode), c s,j is the lithium concentration inside the particle, K s,j denotes subdiffusion coefficients, r ∈ [0, r j ] is the radial coordinate, and r j is the radius of spherical particle. The orders α j ∈ (0, 1] of fractional derivatives are equal to the dispersion parameters of ion subdiffusion in solid electrode particles. The process described by equation (3) implies zero gradients of concentration inside electrodes in the past. It means that the process starts after long relaxation when the distribution of lithium in the electrode particles became uniform, c s,j (t < 0) = c 0 s,j . The boundary conditions have the form: The first condition implies the absence of current in the particle center, and the second one assumes that the current through the electrode surface is equally distributed between the electrode particles. The molar current density is determined by the expression J Li j = I/FA j , where A j is the area of the electrode surface, F is the Faraday constant, and I is the total current flowing through LIB. To write these conditions, the fractional generalization [35] of Fick's law was used.
The current density of ions J Li j is related to overvoltage by the standard Butler-Volmer relation (see, e.g., [36,37]), where the factor J 0,j is defined by the formula J 0,j = k j c max s,j c 0.5 e , parameter k j is a reaction rate constant, the transfer coefficients of cathode and anode charge are chosen to be equal to 0.5, R and T are the gas constant and temperature, and c surf s,j = c s,j r j , t is the lithium concentration near the particle surface. Fraction x j determines the state of charge: From the Butler-Volmer equation, overvoltage can be expressed as follows [32]: LIB voltage is equal to the sum of the differences between overvoltages η j , open circuit voltages U j of electrodes, and voltage drop Φ e,j in the electrolyte: To generalize the results for the case of subdiffusive motion, we use the relation between fundamental solutions of normal and subdiffusion equations. This relationship can be represented in the form [11]: where c nd is the solution corresponding to the normal diffusion equation with a unit diffusion coefficient, g α (t) is the one-sided stable Lévy density, and K α is the subdiffusion coefficient. The relation can be derived from the representation of the subdiffusion process in terms of a subordinated Brownian motion directed by a fractional Poisson process [23].
Using the solution obtained within the standard SP model and presented, for example, in [31], where parameter µ j = J Li j R j /(c max s,j K s,j ), and λ k are the roots of algebraic equation sin λ k = λ k cos λ k , we arrive at the expression for the subdiffusive case with characteristic exponent α j : which includes the Mittag-Leffler function playing a central role in the fractional differential calculus, The Mittag-Leffler function in our solution arises as a result of applying Relation (5) to Solution (6) involving the integral representation: The gamma function appears from the expression for the negative moment of the one-sided stable random variable: By analogy with the normal-diffusion case [38], if we limit ourselves to a parabolic profile, the concentration near the surface is approximately equal to: For given fractions x surf j , we determine overvoltage and open-circuit voltage and then the LIB voltage according to (4). Expressions for the open-circuit voltage U j as a function of the state of charge can be obtained by fitting experimental curves. Expressions for the functions U j (x j ) for Li x CoO 2 and graphite with spherical particles under normal conditions are given, for example, in Appendix A (Formulas (A1) and (A2)) of [31].
In the simplest versions of the SP model, the potential distribution inside the electrolyte is not considered. To model a voltage drop over the electrolyte, the ohmic resistance R cell is chosen: where I is the current density in LIB. The work in [31] presented the approximation of R cell as an empirical function of the ambient temperature and current in LIB. In [39], resistance R cell was modeled as a function of the ionic conductivity of the electrolyte and the thicknesses of the electrodes. At high charge or discharge rates, the potential distribution inside the electrolyte cannot be neglected [40].
To take the distribution of ion concentration across the thickness of the electrodes into account, it is necessary to solve the transport equations in the electrolyte. Under the assumption about ion subdiffusion in the x direction due to the percolative nature of diffusion, we write: We assume normal diffusion in the separator, ∂c e,s ∂t Here, K eff e,j = K e ε brug j are effective subdiffusion coefficients for ions in electrolyte, ε j are the porosities of the anode, separator, and cathode, and brug denotes the Bruggeman coefficient.
The boundary conditions at the separator-electrode interfaces correspond to continuity of ion concentration and current. At the cell boundaries x = 0 and x = L, the "reflecting" boundary conditions are applied [32].
Solutions of these equations can be found in the form of polynomials. The orders of these polynomials for electrodes and a separator may vary. Often, the distribution of electrolyte concentration is approximated by a parabolic profile in electrodes and a linear profile in a separator. We use Formulas (14)-(16) from [32] and modify them using the above procedure using Expressions (5), (8), and (9): The expressions for the coefficients a i and b i are given in [32] and depend on the thickness and porosity of the cathode, anode, and separator and the ion diffusion coefficient in the electrolyte. Due to the cumbersome nature of these expressions, they are not given here. According to the work [32], we calculate Φ e,p (t) x=L − Φ e,n (t)| x=0 : , κ i is electrolyte conductivity in the cathode, anode, and separator. After substituting Solution (10) into the expression for the voltage (4), it is possible to calculate the galvanostatic charge-discharge curves in the subdiffusion version of the SP model. Figures 3 and 4 present the calculation results for different values of the subdiffusion parameters. The basic model parameters are listed in Table 1. The parameters were taken from [16,41] for LIB with the LiMn 2 O 4 spinel cathode, graphite-based anode, and electrolyte LiPF 6 in 1:2 ethylene carbonate and dimethyl carbonate (EC:DMC). The discharge current density was chosen equal to −15 A/m 2 . Decreasing the value of the subdiffusion indices α p and/or α n led to a faster LIB discharge. The new aspect of the degradation of battery performance may be associated, among other factors, with a change in the type of diffusion in LIB components. To verify this statement, let us consider the effect of anomalous diffusion on the LIB impedance and discuss the possibility of determining the anomalous diffusion parameters in individual components using electrochemical impedance spectroscopy.  Table 1 and different values of the subdiffusion index in a cathode.   Table 1 and different values of the subdiffusion index in an anode.

Cell Impedance in the Subdiffusion Model
One of the electrochemical models describing LIB impedance spectra was proposed by Abraham, Kawauchi, and Dees (AKD) in [42]. The model exploits the one-dimensional reaction-diffusion equations, and diffusion is assumed to be normal. The general methodology of the AKD model follows the Newmann approach [43], where volume-averaged transport equations for concentrated solutions are used to describe ion transfer in porous separator and composite electrodes. Using expressions for impedances derived by Abraham et al. [42], one could follow changes in impedance spectra due to variations of physical parameters such as electronic conductivities of electrodes, lithium diffusion coefficients, porosities, etc. The values of the basic parameters used in this work for the AKD model are listed in Table 2. The parameters were taken from [16] for LIB with the LiCoO 2 cathode, graphite-based anode, and electrolyte LiPF6 in 1:1 ethylene carbonate and diethyl carbonate (EC:DEC). The values of the subdiffusion coefficients were taken close to their normal diffusion analogues. The effects of parameter variations on Nyquist plots are presented in Figure 5. Noteworthy is that the apparent slope of the low frequency (quasi-rectilinear) part at the Nyquist diagrams changes with variations of the lithium diffusion coefficient and porosity, which can be interpreted as a sign of anomalous diffusion [6]. That Nyquist plots at low frequencies have slopes different from 45 • could be explained by the joint contribution of ion diffusion in different components (cathode and anode particles, separator, electrolyte, SEI) characterized by different coefficients. As a result, we observed transitions between different normal diffusive processes prevalent in different frequency ranges. The case of "non-ideal" diffusion behavior, when straight sections on Nyquist diagrams have a slope that differs from 45 • , was successfully explained in [6,44] within the framework of the impedance theory of anomalous diffusion [7]. For the sake of fairness, it is worth noting that the low-frequency parts of Nyquist plots observed by [6] were straighter than shown in Figure 5. Therefore, the anomalous diffusion interpretation is a plausible hypothesis. Here, we want to coordinate fractional differential generalization of the SP model with an equivalent circuit model and related impedance spectra.
The Fourier transform of Equation (2) leads to the expression: In this case, the modification of the impedance in the case of subdiffusion is not limited to a simple replacement (iω) on (iω) α : the standard Fick law in the case of anomalous diffusion is not satisfied. For the diffusion Equation (2), the expression for the flux density is the following: In particular, the subdiffusion generalization of Warburg's impedance for a semi-infinite medium has the form [7]: For spherical particles, the Fourier transform of the subdiffusion Equation (3) leads to the differential equation: (iω) α jc s,j (r, ω) = K s,j 1 r 2 ∂ ∂r r 2 ∂c s,j (r, ω) ∂r , 0 < α j ≤ 1, 0 < r < r j .
Using the standard procedure (see for example [32]) from Equation (14) under the relevant boundary conditions and applying the generalized Fick law, one can obtain a generalization of Faraday's reaction impedance: It has the form of an impedance corresponding to a series connections of charge transfer resistance R ct,j and impedance associated with subdiffusion in the solid phase,  Assume ion subdiffusion with exponents β j in electrolyte filling the j electrode (cathode or anode). The impedance of a double layer formed on the surface of a spherical particle can be found by linearizing the transport equation using the Fourier transform: (iω) β jc e,j (r, ω) = K e,j 1 r 2 ∂ ∂r r 2 ∂c e,j (r, ω) ∂r , r > r j .
Solving this equation, in analogy with [45], we arrive at the subdiffusion impedance of a double layer on a spherical particle: The latter expression can be considered as an impedance of parallel connection of a constant phase element and a resistor: The impedance of the anode film is found as the subdiffusion impedance of the spherical layer according to the procedure presented in [45], but using the generalized Fick law. The result represents a parallel connection of the film resistance R f and the subdiffusive generalization of Warburg open impedance: where γ and K f are the dispersion parameter and subdiffusion coefficient for the anode film. Here, however, it should be noted that SEI thickness is usually several nanometers, and the diffusion limit may not be fulfilled on such small scales, while the validity of Expression (16) remains in question.
Assuming independent parallel flow of the charge current due to an electric double layer formation and intercalation current on the active material surface according to [46], after arranging the impedances of the electrodes, separator, and current collector, we arrive at the equivalent circuit shown in Figure 6. This circuit summarizes the popular LIB models (see, for example, [3,32]). Under certain simplifications, the equivalent circuit (EC) presented in Figure 6 can be reduced to the model considered in [3]. Our scheme is also consistent with the one proposed in [32]. However, Li et al. [32] introduced "by-hand" three CPE, substantiating them with the roughness of electrode particles and the fragmented SEI structure. Without denying the plausibility of these arguments, we only note that in our case, the circuit was fully derived within a methodologically-consistent model assuming subdiffusion in the electrolyte filling porous medium and in electrode particles. Consider the simplified circuit presented in Figure 7. The Nyquist diagrams are also presented for four pairs of subdiffusion index values (α n , α p ): (1, 1), (1, 0.5), (0.5, 0.5), and (0.5, 1). For the other parameters: L = 0.12 µH, R I = 0.31 Ohm, C DL,p = 0.026 F, C DL,n = 0.0066 F, R 0 = 0.12 Ohm, R DL,p + R ct,p = 0.11 Ohm, R DL,n + R ct,n = 0.026 Ohm. It can be seen from these diagrams that ion subdiffusions in the cathode and anode have a joint effect on the slope of the low-frequency section. Thus, the slope of this section cannot uniquely determine the subdiffusion indices α n and α p separately from each other. Adverse reactions and degradation processes can lead to a number of undesirable effects associated with the loss of LIB capacity. Usually, LIB aging occurs due to a variety of processes and reactions that simultaneously occur in different areas of LIB. A decrease in LIB capacity is often associated with processes occurring in a negative graphite electrode, where the formation of SEI film leads to an irreversible loss of cycled lithium, to an increase in internal resistance, and to a decrease in the pore volume available for the electrolyte. In our opinion, another additional factor is possible. It is associated with an increase in the fraction of particles, in which anomalous ion diffusion is predominant. The calculated discharge curves in galvanostatic mode (Figures 3 and 4) indicate the plausibility of the hypothesis. An increase in the role of anomalous diffusion due to the disordering of the crystal lattice or the formation of an amorphous film can be estimated by introducing parallel-connected elements (see Figure 8, inset). Suppose that during the cycling, the ε fraction of particles characterized by anomalous diffusion of lithium in the solid phase increases. Let us choose the cathode as an example. The Nyquist diagrams for different values of ε are shown in Figure 8. The model parameters are the same as for Figure 7 with R 1 DL,p + R 1 ct,n = 0.024 Ohm. The diffusion coefficients are taken from Table 1. Elements Z α p s,p and R DL,p are multiplied by ε −1 and the elements Z 1 s,p and R 1 DL,p by (1 − ε) −1 . It can be seen from these diagrams that the slope of the rectilinear segment at low frequencies does not always unambiguously determine the subdiffusion index α and can be both larger and smaller than the slope corresponding to normal diffusion.  Figure 8. Nyquist diagrams for different values of the ε fraction of cathode particles characterized by ion subdiffusive transport. The 1 − ε fraction of particles is characterized by normal diffusion.

Conclusions
In LIB, there are many factors that can sustain the anomalous-diffusive behavior of ions, such as the percolative structure of porous electrodes, an amorphous solid-electrolyte interface, grain boundaries in electrodes and the electrolyte, etc. In our work, the influence of subdiffusion on the galvanostatic discharge and LIB impedance was evaluated within the generalization of the single-particle model. This generalization was based on the fractional differential equations of ion diffusion in electrode particles, SEI, and electrolyte. The distribution of lithium ions in electrolyte and electrode particles was expressed via the Mittag-Leffler function and the Lévy stable density. It was shown that subdiffusion in the solid phase of electrodes and in electrolyte deteriorated the performance of LIB. The calculated discharge curves indicated an acceleration of voltage drop at lower values of subdiffusion index. Within the framework of the new model, an equivalent circuit was derived. It generalized the popular models of LIB. The slope of the rectilinear part of the Nyquist diagram at low frequencies did not always unambiguously determine the subdiffusion index α and could be either larger or smaller than the slope corresponding to normal diffusion. It was shown that the degradation of battery properties can be associated, among other factors, with a change in the type of diffusion in LIB components.

Acknowledgments:
We are grateful to our colleagues Sergey Novikov, Andrey Somov, Pavel L'vov, Ekaterina Morozova, Pavel Kapustin, and Dmitry Vostretzov from the Ulyanovsk State University for the discussion of particular problems related to this work.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: