Thermodynamically consistent derivation of chemical potential of a battery solid particle from the regular solution theory applied to LiFePO4

The chemical potential of lithium in LixFePO4 active cathode nanoparticles and the surface free energy between LixFePO4 and electrolyte were determined with the novel thermodynamically consistent application of the regular solution theory. Innovative consideration of crystal anisotropy accounts for the consistent determination of the dependency of the chemical potential on the mechanistically derived enthalpy of mixing and the phase boundary gradient penalty. This enabled the analytic, thermodynamically consistent determination of the phase boundary thickness between LiFePO4 and FePO4, which is in good agreement with experimental observations. The obtained explicit functional dependency of the surface free energy on the lithium concentration enables adequate simulation of the initiation of the phase transition from FePO4 to LiFePO4 at the surface of active cathode particles. To validate the plausibility of the newly developed approaches, lithium intercalation into the LixFePO4 nanoparticles from electrolyte was modeled by solving the Cahn-Hilliard equation in a quasi-two-dimensional domain.

A rapid increase in mobile applications has resulted in intense research on energy storage materials. Due to their large capacity, Li-ion batteries with several possible active cathode materials are very popular. Maximizing the energy and power densities of such materials while ensuring a long cycle life inherently calls for an enhanced understanding of underlying physicochemical processes. This objective recently resulted in intensified research in the area of modeling of lithium intercalation into active cathode material from electrolyte. Phase field models of Li intercalation are based on solving the continuity equation with the chemical potential derived from the regular solution theory. Due to its phase separation nature [1][2][3][4][5] , Li x FePO 4 (LFP) material is among the most challenging, and it is also one of the most extensively modeled cathode materials. The newly developed theory proposed in this paper was therefore derived for and demonstrated on the phase separating LFP material.
Selected state-of-the-art models 2,4,6-15 offer profound insights into the intercalation process of the LFP, whereas these models feature certain deficiencies related to consideration of the crystal lattice anisotropy and description of the surface-related phenomena at the interface of the cathode material and the electrolyte. Models of the intercalation of lithium into the active cathode material are based on the regular solution theory 16 , which considers intercalated lithium atoms and vacant lithium intercalation sites in the crystal lattice of the cathode material as two diffusing species in the solid solution [17][18][19] . Application of the regular solution theory to such a system is performed by the attributing the solute function to lithium atoms and vacancies and considering the crystal lattice of the cathode material as solvent 16,18,19 . Originally, regular solution theory was designed for the description of liquid solutions with amorphous material structure 1,16,[20][21][22] . This theory assumes an isotropic nature of the considered media, which is not the case for the crystalline cathode material. The obtained chemical potential is thus independent of the crystallographic properties. This deficiency is usually dealt with in such a way that anisotropy of the crystal lattice is introduced into the model by empirical assignation of a tensor nature to the diffusion coefficient and governing parameters of the chemical potential 6,8 . Even though the approach proposed in 6,8 allows for modeling of the lithium intercalation into the active cathode material, the number of independent parameters of the model can be further decreased and thermodynamic consistency of the model can be further enhanced if anisotropy of the crystal lattice is introduced into the regular solution theory.
The second challenge that arises from the use of the regular solution theory for modeling Li intercalation into the active cathode material is related to the consideration of the active particle surface. Phase field modeling of Li intercalation dynamics unavoidably calls for four different boundary conditions 11 . At the phase boundary between electrolyte and active cathode material, the natural boundary condition is applied that interrelates the surface concentration (c) and the boundary surface free energy (γ) 23 . Two common approaches for implementing natural boundary conditions can be found in the literature. The first approach assumes constant surface energy, γ 6,7,10 , which is the simplest possible solution of the boundary problem. Nevertheless, recent experimental and theoretical studies 23 show that such a description of the phase boundary between the cathode material and electrolyte is too oversimplified to efficiently describe all surface phenomena observed over the entire state of charge (SoC) range of the particle. Initiation of the phase transformation from Li-poor to Li-rich phases (nucleation) in phase separating materials mainly occurs at the particle surface, which cannot be adequately described by constant γ 23 . Cogswell et al. 23 showed that with the use of non-constant surface energy, γ, in the model, the nucleation of phases on the surface can be described. This is in line with the second approach of natural boundary condition implementation that assumes linear dependency of γ on c over the entire SoC range of the particle. Non-constant γ leads to initiation of a spinodal decomposition reaction on the active particle surface, which corresponds to the experimental observation 24 . Nevertheless, the linear dependency of γ (a valid description of salt solution surfaces 25,26 ) does not offer a fully thermodynamically consistent description of the interaction between the bulk and the surface of active cathode material 27 . This is because governing equations of the bulk dynamic are, unlike those for the surface dynamics, derived from the regular solution theory.
As regular solution theory plausibly describes lithium dynamics in the particle bulk [17][18][19] , elaboration of more adequate functional relations of γ represents an important contribution towards higher fidelity battery models. Nadkarni et al. 28 (uploaded to arXiv; currently in recension) were the first to tackle this challenge in LFP material by the use of higher order polynomial dependency of γ vs. c. Functional dependency was obtained empirically by polynomial interpolation of the surface free energy values predicted by DFT calculations 29 . Therefore, the proposed methodology relies on a relatively complex workflow.
To provide a solution to the outlined challenges, this paper offers an innovative derivation of the chemical potential of Li in the active cathode nanoparticle from the regular solution theory. This derivation enables a mechanistic consideration of the anisotropy of the crystal lattice in parameters of the governing equation for the chemical potential. Moreover, application of the same theory enables derivation of thermodynamically consistent functional dependency of γ vs. c. Elaborated equations for the chemical potential and the surface free energy can be directly integrated into a particle model, which ensures a high level of predictiveness of the model. To validate plausibility of the newly developed approach, the derived theory was implemented for the LFP cathode material. The original contributions of the paper are therefore as follows: • Explicit derivation of the chemical potential of Li in the active cathode nanoparticle and the interaction energies between diffusing species based on the regular solution theory by considering crystallographic parameters of the cathode material, which also enables mechanistic determination of the gradient penalty coefficient from the mixing enthalpy. • Mechanistic derivation of phase boundary thickness between Li-rich and Li-poor phases in a phase separating cathode material. • Thermodynamically consistent derivation of the surface free energy at the interface between the surface and the bulk of the particle from regular solution theory, where the surface of the active cathode particle is treated as a surface of the regular solution. • A novel quasi-two-dimensional LFP particle modeling framework for modeling Li intercalation dynamics with the boundary condition and chemical potential derived from the regular solution theory.

Methods
Governing equations. Cahn-Hilliard Equation. The temporal and spatial evolution of Li concentration inside active cathode particles is described by the Cahn-Hilliard equation. In general, the Cahn-Hilliard equation can be written as 1,10,30 : where k b T is the thermal energy; ∇ is the vector operator of spatial derivatives; c represents the non-dimensional concentration of lithium (usually non-dimensional molarity), which depends on two independent variables [time (t) and position (r)]; D is the diffusion coefficient, which is generally a function of c; and μ represents the chemical potential of lithium. The chemical potential μ inherently depends on c. In the limiting case of a constant diffusion coefficient and logarithmic dependency of the chemical potential on the concentration, the Cahn-Hilliard equation simplifies to Fick's second law of diffusion. In the case of solid crystalline cathode material, the chemical potential is derived from the regular solution theory. The most general form of the chemical potential reads 2,30 : Scientific RepoRts | (2019) 9:2123 | https://doi.org/10.1038/s41598-019-38635-2 where Ω is the enthalpy of mixing of diffusing species, κ is the gradient penalty between two phases in phase separation materials, and B is the phase boundary strain. The regular solution enthalpy of mixing, Ω, and the gradient penalty parameter, κ, can be obtained from the regular solution theory (e.g. 1,2,[17][18][19][20], and the strain coefficient B can be obtained from the theory of elasticity (e.g. 10,30 ). However, the present derivations of parameters Ω and κ did not take into account properties of the crystal lattice. In these approaches, crystallographic properties of the cathode material are introduced to the model empirically through the use of an anisotropic diffusion constant and the strain term B. This deficiency of present models is approached in the next section, where a thermodynamically consistent derivation of Ω and κ from the regular solution theory that considers crystallographic parameters of the cathode material is presented.
Enthalpy of Mixing and Gradient Penalty. The Cahn-Hilliard equation (Eq. 1) and the chemical potential (Eq. 2) are obtained from minimization of the total free energy functional, which can be derived from the regular solution theory (Appendix). The theory of non-uniform regular solution derived by Cahn and Hilliard 20 can be applied to the crystalline cathode material by considering lithium atoms and vacancies in the crystal lattice as two diffusing species in the solid solution 19 . Regardless of the high anisotropy of cathode material crystal lattices, regular solution (designed for isotropic media in 16 ) allows for description of lithium intercalation dynamics, which is also in agreement with the general consensus communicated in influential refs [17][18][19] , indicating that the regular solution description of anisotropic crystalline material is sufficiently relevant for elaborating state-of-the-art Li intercalation models.
Equations for the chemical potential and concentration dynamics of non-uniform regular solutions (Eqs 1 and 2) were first derived by Cahn and Hilliard 20 , and they led to the determination of parameters Ω and κ. Applied to crystalline cathode material, the original expressions of Ω and κ read [17][18][19] : where c m denotes the maximal molarity of lithium in the cathode material, N represents the areal density of Li intercalation sites, z is the number of intercalation sites nearest to the reference site, and N A is the Avogadro constant. Parameters ε LiLi , ε Liv and ε vv denote pair interaction energies of Li intercalation sites occupied by two lithium atoms, a lithium atom vacancy, and two lithium atom vacancies, respectively. In the theory of Cahn and Hilliard 20 (Eqs 3 and 4), which was developed for isotropic solutions for application on liquid binary solutions, the quasi-crystal approximation was used 16 . The quasi-crystal approximation was first proposed by Guggenheim in the original ref. 16 that postulated regular solution theory; since then, this approximation has been quickly adopted for describing amorphous media. An amorphous isotropic liquid regular solution is treated as an ordered closely packed crystal lattice with a well-defined number of nearest neighbors z and a crystal plane areal site density N. Two parameters of the crystal lattice are introduced to the chemical potential (N and z), although the theory is designed for non-crystalline media 20 . An additional parameter, λ, is needed for the mathematical consistency of the theory 20 . Parameter λ denotes the distance used for defining the concentration gradient in discretized media of the crystal lattice (usually a few inter-particle distances 20 ). In the phase separating materials (i.e., LFP), λ can be interpreted as the thickness of the phase boundary between Li-rich and Li-poor phases 4,10,11,20,30 . From Eqs 3 and 4, the relation between Ω and κ is obtained: m a 2 which is also used by the authors of refs 11,30 . In multiple refs 4,10 , parameter κ is determined by fitting parameter λ to the experimental observations.
As exposed in the Introduction, an innovative approach to deriving the chemical potential from regular solution theory includes explicit thermodynamically consistent determination of parameters Ω and κ, based on crystallographic properties and the interaction energies between diffusing species (lithium atoms and vacancies). This mechanistic evaluation of parameters Ω and κ from the partition function is presented in the Appendix: Eqs A.6 and A.7.
For parameter Ω, the same expression was obtained as that in ref. 20 (Eq. 3), and parameter Ω thus explicitly depends on two crystallographic parameters: N and z. Unlike in other references (e.g. 4,10,11,30 ), an innovative expression is obtained for parameter κ: This expression includes parameter d, which represents the distance between crystal planes, and a newly introduced parameter m. Parameter m is explained in detail in the subsection Evaluation of Parameter m; similarly, to parameters N and d, it depends on the geometry of the crystal plane family under consideration. The proposed approach therefore adequately reflects the crystal structure, including its potential anisotropy, in the parameters Ω and κ, and via Eq. 2, the crystallographic parameters are also adequately reflected in the chemical potential, which is an innovative contribution of this paper.
Resultingly, by combining Eqs 3, 5 and 6, parameter λ can be further determined as: Parameter λ depends on two parameters d and m (see subsection Evaluation of Parameter m), which are mechanistically derived from the crystal properties of the cathode material.
Boundary Conditions. The system of two equations, Eqs 1 and 2, represents a fourth-order partial differential equation system, which unavoidably calls for four boundary conditions: surface where F denotes the lithium flux, subscript bulk refers to the geometrical center of a particle, subscript surface refers to the particle surface, and F BV denotes the Butler-Volmer flux, which is described by the Butler-Volmer equation 31,32 , the standard phenomenological description of electrode kinetics [33][34][35][36] . Two boundary conditions in the particle bulk (Eqs 8 and 9) are obtained from the symmetry of the problem. They are written for the geometrical center of the particle, where reflection symmetry must be satisfied in all dimensions. The first boundary condition for the particle surface (Eq. 10) consistently describes lithium flux through the particle surface from electrolyte. The second boundary condition on the particle surface (Eq. 11) includes surface effects related to the surface free energy of the particle surface in contact with the electrolyte. It is called the natural boundary condition. The challenge in determining this boundary condition (Eq. 11) arises from the unknown functional dependency of the surface free energy between particle and electrolyte on the lithium concentration. One of the innovative contributions of this paper thus arises from the explicit thermodynamically consistent derivation of the natural boundary condition in Eq. 11, which allows for its direct integration into the particle model.
Natural Boundary Condition. The natural boundary condition on the surface between LFP and electrolyte is obtained by minimization of the total free energy functional with application of the Lagrange variational principle 37 (Appendix). The natural boundary condition implies proportionality between the derivative of the surface-energy concentration (∂γ/∂c) and the concentration gradient ∇c on the boundary (Eq. 12). The proportionality coefficient is equal to the magnitude of the gradient penalty term κ: where ∂ denotes the boundary of the calculation domain (surface of the particle). Two simpler approaches with the constant γ value and the linear γ(c) dependency, which were addressed in the Introduction, result in the following boundary conditions 11 : where β is a constant.
Constant γ ∂ ∂c and Hildebrand formula. As addressed in the Introduction, the stability and accuracy of the models can be improved if the particle surface is considered as the surface of a regular solution. The surface of the regular solution was extensively studied in the case of liquid regular solutions by several authors [38][39][40][41] . Hildebrand showed from the Gibbs adsorption formula that the surface energy of the regular solution must satisfy the following condition 27 : where c * denotes the equilibrium bulk concentration and Γ is adsorption, which measures how much lithium is adsorbed on the surface. The Hildebrand relation (Eq. 14) is not satisfied for the most commonly used linear dependency of γ(c) (used in the boundary condition in Eq. 13), which reads: FePO LiFePO 4 4 where γ denotes the surface free energy of the mixture, γ FePO 4 and γ LiFePO 4 represent free surface energies of pure phases FePO 4 and LiFePO 4 , and c denotes the lithium surface concentration. This leads to a conclusion that a more complex model of surface concentration dynamics should be incorporated in order to provide full thermodynamic consistency. The special case of Eq. 15, where γ γ = FePO LiFePO 4 4 , results in constant γ(c), and consequently, ∂γ/∂c = 0. This γ(c) dependency satisfies Hildebrand's condition for a regular solution in the case where no adsorption on the surface occurs, whereas it is subjected to the deficiencies addressed in the Introduction.
Surface free energy of the regular solution. The functional dependency of the surface free energy γ(c) of an active cathode particle surface in contact with electrolyte was derived by the Murakami approach 40 . Murakami et al. 40 proposed derivation of the regular solution surface energy as a generalization of better known approaches of Guggenheim 38 and Defay and Prigogine 39 . It is derived for the general case of an infinite surface of isotropic regular solution of the solvent with two mixing solutes and a vacuum. Murakami et al. 40 used the quasi-crystal approximation for the liquid solution description, which approximates an amorphous liquid phase as an ordered quasi-crystal with well-defined positions of nearest neighboring atoms of the reference atom. Therefore, their theory can be generalized and applied also to the ordered crystalline solid solution.
An original contribution of this paper, addressed in the third bullet point of Introduction arises from the first successful application of the Murakami et al. 40 approach, which was elaborated for the isotropic media that takes in consideration only nearest neighboring atoms, to the anisotropic crystal lattice. To consistently account for the phenomena in batteries, governing equations were derived for the crystal solid solution of lithium in LFP in contact with electrolyte. This resulted in explicit dependency of the surface free energy of the phase boundary between LFP and the electrolyte on the lithium surface concentration, given in the Appendix (Eq. A.5), since the areal integral part of the total free energy F is by definition equal to the surface free energy density γ. The functional dependency of the surface free energy between LFP and electrolyte thus reads: The first two terms in Eq. 16 represent the linear dependency, which is the same as in Eq. 15. In addition to the first two linear terms, Eq. 16 includes a quadratic Ω term, multiplied by parameter m (defined and explained in the next section). This term in Eq. 16 can be interpreted as a correction to the mixing enthalpy of the surface layer due to absent nearest neighbors. Higher-order terms in Eq. 16 arise from the changed ratio between the type of interacting species pairs (i.e., pair ratios of Li-Li, vacancy-Li and vacancy-vacancy) due to adsorption of lithium to the particle surface and the interaction of diffusing species with electrolyte. Figure 1 represents the plot of surface free energy as a function of surface concentration, where separate contributions of different terms of Eq. 16 are presented. Parameters used in the Fig. 1  is taken from the reference of Ferguson et al. 23 . Figure is   The red line in Fig. 1 represent the first two terms of Eq. 16, which represent the linear dependency of γ(c) (Eq. 15). The green line in Fig. 1 shows the first two linear terms and the quadratic contribution of the third term. The blue line in Fig. 1 shows a full dependency of the surface free energy on the concentration in Eq. 16. The differences between the linear γ(c) function (red line in Fig. 1) and the newly derived γ(c) dependency (blue line in Fig. 1) is most pronounced at high and low surface concentration (beneath the first spinodal point and above the second spinodal point). Therefore, the new thermodynamically consistent derived γ(c) dependency most severely enhances the lithium dynamics near the surface in the solid solution regime, since the slope of the function γ(c) is most severely altered by the correction terms in these two regions. This enables consistent description of Li surface dynamics over the entire SoC range, since newly derived explicit γ(c) dependency does not initiate the phase transition in LFP prior to reaching the spinodal concentration, which is in agreement with experimental observations 23,24 .
Evaluation of Parameter m. Equations 5, 6 and 16 include the explicit dependency on a newly introduced parameter m (derivation in Appendix). Model parameter m was determined from the crystallographic properties of LiFePO 4 .
In the classical quasi-crystal formulation of regular solution theory, only interactions between nearest neighbors are taken into account [38][39][40] . When applying regular solution theory to the description of the solution surface, absent nearest neighbors on the surface are the main contribution to the surface energy 38 . The regular solution model (Eqs A.1 and A.3 in Appendix) accounts for this contribution through the mechanistically based parameter m. m represents the half fraction of the nearest neighbors outside the crystallographic plane parallel to the surface, and its product with the number of nearest neighbors coincides with the number of absent nearest neighbors at the surface. m is thus a measure of the deficit in interaction energy due to the absent nearest neighbors at the surface.
Due to the high anisotropy of the LiFePO 4 crystal lattice (Fig. 2), every Li site in the crystal only has two nearest neighbors (in b direction, Fig. 2). The contributions of all other Li crystal sites to the interaction energy are thus not taken into account if classical quasi-crystal formulation of regular solution theory is applied to the anisotropic crystal lattice.
To solve this challenge, a modification of the existing approach is proposed in this paper. The meaning of parameter m from the quasi-crystal formulation of regular solution was generalized in order to consider contributions of all Li crystal sites to the total interaction energy. To preserve the main significance of parameter m (i.e., a measure of the deficit in interaction energy due to the absent nearest neighbor at the surface), an innovative generalization of m for the anisotropic lattice was designed. Parameter m was determined as a sum of normalized individual contributions of all neighbors that lie outside of the lattice plane parallel to the surface (lattice sites marked blue in Fig. 3(b)).
Li intercalation sites (yellow spheres in Fig. 3(a)) themselves form a hexagonal crystal lattice. Pair interactions that depend only on the distance between two particles can be calculated for sites of such lattices. The magnitude of the contribution of an individual Li site to the interaction energy decreases with distance. The total interaction energy is obtained by spatial integration of the pair interaction potential (φ(r)) multiplied by the pair distribution function (w(r)) 42  The calculation of normalized contributions of each Li-intercalation site to the total interaction energy shows that a non-negligible contribution to the pair interaction comes from pairs that are less than 0.61 nm apart. The contribution of the pairs that are 0.61 nm apart is on the order of a percent, whereas the contribution of the next furthest pairs is an order of magnitude smaller. Every Li site has fourteen neighbors within a 0.61 nm radius ( Fig. 3(a)), which is the upper summation limit in Eq. 18. The normalized contributions of all fourteen nearest neighbors to the interaction energy are listed in Table 1.
m was calculated for the surface parallel to the crystal planes with the orientation of (101). This is because (101) is the crystal orientation that shows the lowest misfit between LiFePO 4 and FePO 4 crystals 10 (lowest parameter B in Eq. 20). As a consequence, all phase boundaries between LiFePO 4 and FePO 4 are oriented normal to the 101 direction due to strain energy minimization 10 . From the projections in Fig. 3(b), the type and number of Li-site neighbors inside and outside the crystal plane parallel to the surface can be determined, which defines parameter m (blue and green spheres in Fig. 3(b)) as: 19 is the number of atoms that lay outside the crystal plane for each of five crystallographic direction of neighbors inside 0.61 nm distance (determined from the Fig. 3). For the calculated (101) orientation of the LFP that is simulated in proposed single particle model (Eqs 20, 21 and 23), m = 0.101. The elaborated methodology is applicable also to other crystallographic orientations of the surface. Table 2 shows values of parameters m for five differently oriented surfaces.

Modeling
To validate the plausibility of the newly developed approach for determining the chemical potential based on the regular solution theory by considering the crystallographic parameters of the LFP cathode material (Eqs 2, 3 and 6) and the functional dependency of surface free energy on lithium concentration (Eq. 16), a simulation model was developed. Phase field modeling of lithium intercalation was performed in one dimension ((101) crystal direction in LFP) with additional discretization of surface planes in the perpendicular direction ((010) crystallographic direction) for accurate implementation of newly derived surface phenomena (the discretization of the  computational domain is schematically shown in Fig. 4). This forms a basis of a quasi-two-dimensional model that represents a good trade-off between model complexity, computational demand and relevance of the results of the model. Nevertheless, the derived methodology is universal, and it can be used for modeling Li intercalation into LFP within a modeling framework of arbitrary dimensionality. The Cahn-Hilliard equation was solved in the two-dimensional domain for the Li intercalation simulation: As elaborated in the Governing Equations section, four boundary conditions were needed for the solution of a fourth-order PDE system. Concentration and chemical potential gradients were set to zero at the center of the particle due to the symmetry of the problem. At the surface parallel to the (101) crystal planes, Butler-Volmer flux in the x direction was defined the gradient of chemical potential. The concentration gradient on the particle surface was determined from calculated γ(c) dependency (Eq. 16). The surface flux across the (010) oriented surface (direction of fast diffusion) was modeled by the source term (Allen-Cahn model 3-5 ) in Eq. 20. For this specific case, the boundary conditions in (Eqs 8-11) are: All model parameters and constants of the model are listed in Table 3.

Results and Discussion
The results are presented as spatial plots of the non-dimensional lithium molarity (c/c m ) and the one-dimensional chemical potential. The x-axis represents the (101) crystallographic direction. All simulations were performed for constant overpotential η. Results are presented on two different particle sizes (100 nm and 250 nm particles), since this are the representative particle sizes in commercial LFP cathodes 43,44 . Many different physicochemical phenomena can be seen from the results presented in Figs 5, 6, 7 and 8. Figure 5 shows a simulation of near thermodynamic equilibrium of Li intercalation. Figure 5(a) shows the non-dimensional Li molarity, and Fig. 5(b) displays the corresponding chemical potential. Figure 5(c-f) show separate terms from the chemical potential (Eq. 21) at different SoCs of the particle. This corresponds to temporal evolution of concentration and chemical potential inside a 100 nm particle for the case of a near-equilibrium lithium intercalation from electrolyte to particle at a low overpotential of η = 1 mV.
Lithium intercalation to the active particle undergoes a solid solution regime until the first spinodal point is reached (Fig. 6)   occurs at the particle surface, which is in agreement with findings in other references (e.g., Cogswell et al. 23 ). The phase separation nature of the LFP material is well represented in Fig. 5(a). After the nucleation of the Li-rich phase at the particle surface, the phase field is decomposed into two regions: regions with high lithium molarity (LiFePO 4 phase) and low lithium molarity (FePO 4 phase).
The simulation in Fig. 5 was conducted for a very low overpotential, which produces near thermodynamic equilibrium conditions. Very low overpotential η yields a very low flux of lithium from the electrolyte to the particle bulk. Both phenomena are clearly reflected also in a very small spatial variation in the chemical potential inside the particle at a particular SoC value ( Fig. 5(b)). Consequently, the temporal evolution of the Li concentration field in a particle follows a trend characteristic for near thermodynamic equilibrium conditions of the LFP material, which for a 100 nm particle (Fig. 5) results in domino cascade dynamics 4,9 , where two phase boundaries travel from the particle surface to the inside of the particle i.e., perpendicular to the fast diffusion direction ((010) crystallographic direction).
To clearly demonstrate spinodal decomposition in a particle, the particle-averaged value of the chemical potential obtained by integration over the entire computational domain is plotted as a function of a SoC of the particle (Fig. 6). The first fact that can be derived from the spinodal curve in Fig. 6(a) is that the phase transformation initiates at spinodal composition, which can be seen as a shallow local minimum at the top of the left spinodal potential barrier (the part of the plot in Fig. 6(a) at SoC = 0.128 marked with red circle). At the moment of LiFePO 4 phase nucleation on the surface, the transient phenomenon of quick lithium adsorption to the surface occurs. A phase boundary emerges inside the particle bulk at this moment, which results in an additional contribution to the chemical potential due to strain and gradient penalty (Eq. 2 and 5(c) and (d)). In a very small SoC window (approximately from SoC of 0.128 to 0.132), a few crystal planes at the (101) surface are far from thermodynamic equilibrium (Fig. 6(a), although the entire particle is in near thermodynamic equilibrium conditions. This observation can be explained also on the basis of a comparison between chemical potential of the bulk and the surface, calculated only from contributions of the entropy and the mixing enthalpy term of Eq. 2. Figure 6(b) shows chemical potentials of the bulk (blue plot in Fig. 6(b)) and the surface (red plot in Fig. 6(b)). Both were calculated from the entropy and the mixing enthalpy, since these two terms are the only ones that are explicitly dependent on c. An additional term that accounts for surface effects was introduced as the surface chemical potential (red plot in Fig. 6(b)) that was derived from Eq. 16. The maximum value of the chemical potential at the surface is reached prior to the bulk (lower value of red curve in comparison to blue curve at SoC = 0.128 in Fig. 6(b)). This explains adsorption of lithium to the particle surface.  At the second spinodal point, collision of intercalation waves 4 can be seen in Fig. 6(a) (marked with a green circle). It can be seen in Fig. 5(a) that at this moment, two intercalation fronts meet at the center of the particle. The gradient penalty and strain terms drop to zero, which leads to a rapid decrease in the chemical potential, which is discernible in Fig. 6.
To investigate the significance of contributions of the gradient penalty, strain, enthalpy of mixing and entropy to the chemical potential (Eq. 21), separate contributions of each term in Eq. 21 are plotted in Fig. 5(c-f)  . Low current Li intercalation in particle with 100 nm length at overpotential η = 1 mV. (a) Nondimensional molarity of lithium spatial dependency, (b) chemical potential spatial dependency, (c) gradient penalty contribution to chemical potential, (d) strain energy contribution to chemical potential, (e) enthalpy of mixing contribution to chemical potential, and (f) entropy contribution to chemical potential. Rescaled parameters used in the axis labels can be written as follows phases was estimated in two different ways. It was calculated for the (101) crystallographic orientation of phase boundaries from Eq. 7. The calculated value was compared to the one obtained from a simulation in which the phase boundary thickness is recognized as the distance between an adjacent minimum and maximum in a gradient penalty contribution to the chemical potential (Fig. 5(c)). A plausible value of λ = 2.7 nm was obtained by both approaches, which confirms adequacy of Eq. 7. The phase boundary thickness thus accounts for approximately five crystal planes. Parameter m can, therefore, also be interpreted as an inverse square of the number of crystal planes in the phase boundary between Li-rich and Li-poor phases (Eq. 7). Figure 7 shows the results for a near equilibrium lithium intercalation at η = 2 mV for a larger 250 nm particle. Similar physicochemical phenomena, as already discernible in Fig. 5, can also be seen in Fig. 7. A phase transition is initiated at the particle surface. After that, Li-rich phase regions advance from the surface to the particle bulk in the domino cascade manner. Due to the larger particle size, in comparison to the particle from Fig. 5, another Li-rich phase region emerges at the center of the particle. In this way, four phase boundaries are created inside the   particle. The newly emerged Li-rich region has the shape of a bend parallel to the (101) crystal plane, which is in agreement with the model of Cogswell et al. 10 . The thickness of Li-rich bands increases with time in the domino cascade manner due to the low overpotential. At a high overpotential of η = 20 mV and a particle size of 100 nm (Fig. 8), the magnitude of the Butler-Volmer term in Eq. 20 exceeds the magnitude of the first diffusion term, and phase separation in the particle during lithium intercalation is suppressed. Li flux from electrolyte to the particle bulk through the (010) crystallographically oriented particle surface is large in comparison to the transverse uphill diffusion (Cahn-Hilliard dynamics in x direction in Fig. 8). Consequently, the phase field of Li concentration remains nearly constant. This coincides well with the results of in-situ X-ray diffraction measurements performed and published by Lim et al. 24 . On the other hand, the spatial dependency of the chemical potential is far from a nearly constant value in this case, indicating the non-equilibrium nature of lithium intercalation at this overpotential. The transition from the solid solution regime to the phase separation regime therefore does not occur in this case.
The presented results thus confirm that the proposed novel thermodynamically consistent derivation of the chemical potential of Li in an active cathode nanoparticle (which includes innovatively derived terms of gradient penalty, phase boundary thickness between Li-rich and Li-poor phases, and surface free energy at the interface between surface and bulk of the particle) allows for adequate simulation representation of recent experimental findings in the LFP material. In addition, the mechanistic basis of governing equations enables the achievement of credible simulation results over the entire SoC range, for arbitrary charging and discharging rates and particle sizes, confirming their general applicability. Furthermore, explicit analytic derivation of all equations allows for their direct integration into the particle model, which is another advantage of the proposed approach. This is also one of the main strengths of the newly developed approach that motivates future applications, as electrodes consisting of phase separating materials are subjected to inter-and intra-particle phase separation being governed by a complex interplay between lithiation levels of particles, particle sizes and their potentials, and thus direct integration of the proposed governing equations in the multi-particle models allows for obtaining new insights also on the electrode level. This is of particular importance for more profound analyses of transport phenomena on the electrode and on the particle level, which are directly related to degradation phenomena driven by mechanical stresses due to inter-and intra-particle phase separation. Such type of applications thus pave the way towards simulation supported engineering of electrodes with enhanced performance and prolonged service life.