Application of phase-field method in rechargeable batteries

Rechargeable batteries have a profound impact on our daily life so that it is urgent to capture the physical and chemical fundamentals affecting the operation and lifetime. The phase-field method is a powerful computational approach to describe and predict the evolution of mesoscale microstructures, which can help to understand the dynamic behavior of the material systems. In this review, we briefly introduce the theoretical framework of the phase-field model and its application in electrochemical systems, summarize the existing phase-field simulations in rechargeable batteries, and provide improvement, development, and problems to be considered of the future phase-field simulation in rechargeable batteries.


INTRODUCTION
Due to the rapid consumption of non-renewable fossil fuels and aggravation of environment problems 1 , energy storage becomes a fundamental issue for the integration of renewable sources into electric power systems 2,3 . Therefore, energy storage plays an essential role in future energy sustainable development 3,4 . Breakthroughs in energy storage technology can make energy distribution and adjustment across time and space, which has revolutionary significance to the production and consumption of energy. During the past decades, the research of electrochemical energy storage has developed rapidly, especially in rechargeable batteries (RBs), has been carried out. Currently, RBs still have many problems, such as high cost and thermal safety problems; hence, it is necessary to improve the energy efficiency, life span, safety, and cost performance, which requires a combination of theory, experiment, and computation to accelerate the research and development of key materials.
To discover, develop, manufacture, and deploy advanced materials at least twice as fast as is possible today, at a fraction of the cost, the US Government announced the Materials Genome Initiative in June 2011 5 . This initiative calls for significant advances in three research areas: computation platform (i.e., multiscale computational materials science), database platform (i.e., opensource cyber infrastructure for data management), and experiment platform (i.e., an integrated approach combining computation and experiments) 6 , aiming to develop more accurate and powerful Integrated Computational Materials Engineering tools to replace expensive and time-consuming experimental procedures. In general, there are three methods in materials sciencetheoretical, experimental, and computational researches. Computational research is neither as difficult as theoretical research nor as time-consuming and laborious as experimental research. Rapid progress of computer technology makes computational research become a significant means for investigating microstructure evolution of material systems. The phase-field model is a computational research model, which describes microstructure evolution of material systems based on thermodynamics. In this review, we mainly focus on application of the phase-field model in RBs, including both the existing simulation results and future development trend.
Theoretical framework of the phase-field model is discussed in 'Basics of phase-field model' section. Existing phase-field simulations and further development and problems for phase-field simulations of RBs are summarized and discussed in 'Existing phase-field simulation results for RBs' section. Concluding remarks and perspectives are given in 'Concluding remarks and perspectives' section.

BASICS OF PHASE-FIELD MODEL
In thermodynamics, a phase is defined as a spatial region with relatively uniform and stable structure, and a field is a function of space and time. Literally, the phase-field model is a computational model which describes microstructure evolution of material systems as a function of space and time. One feature of the phase-field model is the diffuse interface. Slovenian physicist Stefan proposed the sharp interface method around 1890 to deal with phase-boundary problems. The connective conditions between phases are essential for solving the sharp interface model, which further needs tracking of phase boundaries. However, morphology and motion of phase boundaries for practical problems may be rather complex, which brings tremendous difficulty for solving the sharp interface model. The diffuse interface model, i.e., phase-field model, introduces the phase-field variable to indicate the phase state of the system, which varies continuously across a phase boundary. Evolution of the phase-field variable directly reflects the phase transition process, and no phase-boundary tracking is needed, which reduces complexity of the problem enormously. The other feature of the phase-field method is convenience for dealing with applied fields. If the system involves the concentration, temperature, and applied electromagnetic field except phase transition, more phase-field variables are needed to describe these factors of the system. All the field variables introduced form a complete variable set for describing the system, whereby the free energy functional of the system can be constructed. Evolution equations of field variables are derived according to principles of local equilibrium 7 and free energy minimization 8 , which are usually nonlinear partial differential equations. An applied field can be considered as a part contribution to the free energy functional and its effect on evolution of the system can be reflected through its effect on the free energy functional. Microstructure evolution of material systems can be extracted after solving the evolution equations numerically. It can be seen that the description way of phase-field for material systems is consistent with that of Maxwell's theory for electromagnetic fields. In essence, the phase-field model is an application of thought of classical field theory in materials science.
Theoretical framework of phase-field model In the phase-field model, the most basic field variables are phasefield variables ϕ α ðα ¼ 1; ; nÞ, which describe phase states of the system. Phase-field variables may have the different meanings in different phase-field models, but they are always regarded as non-conserved variables since phases can be produced or disappear. Their evolution obeys the Allen-Cahn equation 9 In general, there is not only one component in the system; so concentration variables c i ði ¼ 1; ; mÞ should be introduced to describe the system. If there is no chemical reaction, concentration variables are regarded as conserved variables. Their evolution obeys the Cahn-Hilliard equation 10 In Eqs. (1) and (2), L αβ is the interface mobility between phases α and β, M ij is the atomic mobility between components i and j, f is the free energy density, and δf ∂∇cj are the variational derivatives. Compared with common derivatives, variational derivatives include the contribution of field variable gradients. For solving Eqs. (1) and (2), the total free energy of the system (including the chemical free energy, elastic energy, electric energy, etc.) should be expressed as a functional of the introduced field variables: (3) f elec and f elas can be calculated according to electrodynamics and elasticity theory, respectively. The chemical free energy density of a single phase f chem α can be obtained from computational thermodynamics. However, f chem cannot be obtained directly from computational thermodynamics at present and it can be constructed reasonably through f chem α as follows: where W αβ is the potential barrier, ε α and κ i are the coefficients of additional energy from phase-field and concentration gradients, and h ϕ α ð Þ ¼ ϕ 3 α ð6ϕ 2 α À 15ϕ α þ 10Þ is the weight function of phase α. The phenomenological parameters W αβ , ε α , and L αβ can be evaluated from the relations between them and some physical parameters. It should be pointed that M ij and κ i are not phenomenological parameters of the phase-field model, and they can be determined by thermodynamic and diffusion kinetic databases. Equation (4) has very explicit physical meaning: the chemical free energy density includes the homogeneous part (the former two terms, which are only relevant to field variables) and inhomogeneous part (the latter two terms, which are relevant to field variable gradients). The homogeneous part includes the weighted average of all individual phases (first term of right side) and interface potential barrier (second term of right side), and the inhomogeneous part includes the contributions from phase-field gradients (third term of right side) and concentration gradients (fourth term of right side).
Electrochemical phase-field model Application of the phase-field model in electrochemical systems involves electric field and sometimes involves elasticity. The electric potential Φ and stain tensor ε are (possibly) introduced to describe the electrochemical system. f elec and f elas can be expressed as follows: where ρ elec is the charge density, F 0 is the Faraday constant, z i is the valence of components i, V mol is the molar volume of the system, ε α and ε 0 α are the total and eigen stain tensors of phase α, and C α is the stiffness tensor of phase α. Strictly speaking, elasticity in electrochemistry is usually a kinetic problem. The relaxation process is very quick compared with diffusion and phase transition, so for simplicity, elasticity is treated as an equilibrium problem and it is generally assumed that the stresses or strains of two phases are equal in the interface. The former, known as Reuss limit 11,12 , is applicable to the case that at least one solid phase is soft; the latter, known as Voigt limit 12,13 , is applicable to the case that two solid phases are both hard.
The equations of electric potential and stain are as follows: where σ elec is the conductivity of the electrochemical system. Equations (1) and (2) must be modified for the electrochemical system because of existence of chemical reactions, i.e., the source term n i;α (produced moles of component i in phase α per unit volume and time) should be taken into account. The modified phase-field and concentration equations can be expressed as The source term n i;α may be taken as different expressions in different papers, but it is generally a function of the current density, phase-field, and concentration. The current density is given by the famous Butler-Volmer equation where j 0 is the exchange current density in equilibrium, R is the gas constant, a is the coefficient of activation energy in oxidation reaction, and η ¼ 1 δci is the overpotential. Equations (7)-(10) constitute basic kinetic equations, which describes microstructure evolution of electrochemical systems.  19 developed a continuum phase-field model considering anisotropic (orthorhombic) and inhomogeneous elastic effects to simulate the diffusion-limited phase transformation during Li intercalation in LiFePO 4 nanoparticles. Hong et al. 20 proposed a hybrid mode of phase transformation to explain the phase-growth morphology and kinetics. The results showed that the surface reaction and Li diffusion in different crystallographic directions control the phaseboundary migration 20 . In addition, many researchers have further studied the relationship between Li diffusion and phase transition in LiFePO 4 16,21,22 , phase stability in LiFePO 4 nanoparticles [23][24][25] , and response performance of the LiFePO 4 cathode 26 . For example, Chueh et al. 27 explored the influence of surface diffusion on Li x FePO 4 phase transition. Their results showed that the Li x FePO 4 is a three-dimensional Li-ions conductor with three plane migration paths: bulk diffusion, surface diffusion and electrolyte diffusion 28 .
Solid-state batteries are promising candidates increase the safety and long-term circulation capacity 29 . In 2017, Chen et al 30 . developed a multiscale model which integrates the DFT calculation and phase-field simulation. The effective ion conductivity of β-Li 3 PS 4 electrolyte was predicted via the phase-field simulation (Fig. 1). In 2019, Gamon et al. 31 studied Li-Al-O-S series compounds via the combining the experiment and phase-field simulation. They discovered a sulfide Li 3 AlS 3 , which expands the possibility for exploring new interesting structures via phase-field simulations 31 . In 2020, combining the DFT calculation and phasefield simulation, Shen et al. 32 revealed the role of grain boundary in ionic conduction of solid electrolyte.
Stress evolution and fracture Previous stress-studying phase-field models in the anode usually treated the lithiation stress as a diffusion-induced stress, i.e., Li diffusion in electrodes makes the composition deviates from its stoichiometric state [33][34][35][36][37][38] . If Li distribution is non-uniform, deviation from stoichiometry usually results in a volume change and generates stress [39][40][41] . Prussin 42 dealt with diffusion-induced stress by analogy with thermal stress. In 2005, Garcia et al. 43 first studied the stress due to Li-ion intercalation and they split the total strain into elastic strain and inelastic strain. At the same time, Christensen and Newman 44,45 considered the mechanical stress and fracture during Li-ion intercalation.
Si is a promising next-generation anode material for Li-ion batteries because of the low cost and high theoretical capacity. Lithiation and delithiation in the Si anode result in elastic and plastic deformation 46 . Chen et al. 37 proposed a phase-field model coupled with large elastic-plastic deformation to investigate the phase evolution, morphology, and stress mitigation inside the crystalline Si electrode during Li intercalation. Zhao et al. 35,47 and Bower et al. 33 formulated computational frameworks to model diffusion-induced plastic deformation. Walk et al. 48 combined the Cahn-Hilliard equation with small and large deformations, respectively, to analyze the stress evolution. Gao and Hong 49 further developed a phase-field model, which considered both volumetric and deviatoric inelastic deformation as direct results of the lithiation at the reaction front. Harris et al. 50 studied Li ion transport and Li elastoplasticity at electrode-electrolyte interfaces, and the results showed that contact elastic-plasticity depends on the ratio of Li yield strength to composite roughness [50][51][52] . To reduce the volume change of the Si anode during charging-discharging process, nano-silicon anode materials are mostly used in experiments. In 2014, Yang et al. 53 proposed a chemo-mechanical model to study the deformation of Si nanowire (NW) during lithiation. In 2015, Xie et al. 54 developed a phase-field model to study the structure deformation and stress evolution of Si NW during lithiation. There are large local stresses and stress gradients at the initial lithiation process, indicating that NW structure may generate damage and fracture 54 . Moreover, several models and studies have been reported on the electrochemical and mechanical responses of the amorphous Si electrode during the Li-intercalation process [55][56][57][58] .
Stress generated by diffusion via expansion and shrinkage in intercalation and deintercalation of metal cations causes microcracks, which leads to fracture with large deformation 39,59 . The mechanism of the crack propagation in electrode during charging-discharging process is critical to understand structure stability of the battery electrode 59 71 and Stein et al. 72 conducted a comprehensive study on the effect of chemical reaction on fracture surfaces (Fig. 2). The state of charge is the ratio of amounts of Li ions inside the particle at current state and those at full lithiation state 72 .
Electrodeposition and dendrite growth In Li-ion batteries, the formation of Li dendrite during charging processes is one of the main obstacles in their widespread application, which will result in decrease of reversible capacity and an internal short circuit 28,73 . Upon Li deposition, the surface of the electrode particles usually cracks due to volumetric expansion, and new Li is exposed for further reactions 74,75 . Li dissolution in discharging processes creates pits and crevices with low impedance, and Li-ion flow at the defects leads to rapid growth of metal filaments and dendrites 76,77 . Therefore, understanding the physical mechanisms of this complex non-equilibrium process is essential to improve performance of Li-ion batteries. In 2003, Monroe and Newman 78 first established an electrochemical model for dendrite growth, which is aimed at evolution of dendrite tip height and growth velocity in Li-polymer batteries. The next year, Guyer et al. 79,80 developed a 1D phase-field model to investigate the charge separation on the interface and kinetic behavior of electrodeposition process. Later, Shibuta et al. 81 coupled the Cahn-Hilliard equation with the Butler-Volmer equation to simulate the 2D electrodeposition. Ely et al. 82,83 developed a phase-field model to describe growth kinetics of Li deposits. However, all these models assume a linear kinetics, which are not applicable for systems far from equilibrium. Liang and Chen 70,84 presented a nonlinear phase-field model, which takes the Butler-Volmer reaction kinetics into account, which can simulate and predict the Li deposit growth during charging processes without considering the SEI layer effect. Enrique et al. 85 generalized DeWitt's model to study the morphological evolution of Limetal electrode during electrodeposition 86 . In 2015, Chen et al. 87 established a consistent thermodynamic phase-field model to investigate the dendritic patterns, as shown in Fig. 3.
A plenty of works have been done to research the influence factors of Li dendrite growth. Jana et al. 83 presented a phase-field model to assess the effect of current density on lithium dendrite growth. Yurkiv et al. 88,89 studied the influence of the stress and the SEI on Li electrodeposition, and reproduced filament structure of Li observed in experiments. Harris et al. 90 studied the effect of external pressure on the dendrite growth in lithium metal batteries. Their results showed that if there is sufficient local stress, Li avoids plating at the tips of growing Li dendrite. Hong et al. 91 developed a grand potential-based phase-field model, which reveals the fundamental physical insights related to the intimate, dynamic competition between dendrite growth, ion transport, and electrochemical reactions. Liu et al. 92 discovered that the shape of SEI can limit the dendrite morphology, e.g., denser SEI lead to slower dendrite growth. Recently, Zhang et al. 93 described the Li dendrites growth with various conductive structures by using phase-field model. Xu et al. 94 first studied the adhesion effect of lithium dendrites between the interface layer and lithium metal and they found that enhanced adhesion energy can help inhibit the growth of Li dendrite. Chen et al. 95 investigated the thermal effect on the Li dendrites growth. Hong et al. 96 further revealed that the temperature-dependent performance of batteries is related to the reaction barrier and ionic diffusion barrier.
Besides, Zn is regarded as a promising anode material due to high theoretical capacity, non-toxicity, and abundant sources. The main problems are corrosion and dendrite growth. Cogswell 97 developed a phase-field model based on Marcus kinetics for concentrated solutions to simulate dendritic growth of Zn. In 2015, Wang et al. 98 studied the influence of overpotential, surface energy anisotropy, electric field and ion concentration on the dendrite growth and morphology. In 2019, Wang et al. 99 studied the influence of overpotential, interface energy anisotropy, ion diffusion rate, electrolyte concentration, and electrolyte conductivity on the morphology of Zn electrodeposition. Very recently, Yurkiv et al. 100 investigated ion concentration, electrostatic potential, stress strain, and dendritic morphology of zinc electrodeposition. Hong et al. 101 studied the electrodeposition of Zn and designed a hybrid aqueous ZnSO 4 /solid polymer composite electrolyte to delay the dendrite growth.  Further development and problems to be considered Improvement and development. To our understanding, there are two kinds of electrochemical simulation. One is describing a kinetic phenomenon in electrochemistry, which needs a reasonable free energy functional. The other is a kinetic process of a practical electrochemical system, which needs the realistic thermodynamic database to provide free energy functional of the system. Most existing phase-field simulations of batteries materials belong to the former where their electrolyte solutions are treated as ideal solution. Obviously, such simulations cannot reflect individuality of the investigated systems. Therefore, for phase-field simulations of RBs, future improvement should focus on coupling thermodynamic databases of the electrode and electrolyte to investigate individuality of the systems and its difference with a similar battery system. In thermodynamics, the free energy of single phases is constructed by molar fractions. For multiphase systems, the forgoing molar fractions mean phase molar fractions, which are referred to phase concentrations. However, the diffusion kinetics only provides the evolution equation of total concentrations. Consequently, there are two ways to connect thermodynamic databases and phase-field simulations: (1) using total concentrations to describe evolution of components and decomposing them into phase concentrations by a certain rule, and (2) constructing the evolution equation of phase concentrations. The former is represented by KKS model 102 , wherein the principle of equal chemical potential is used to decompose total concentrations. However, this principle is only applicable for equilibrium and weak nonequilibrium states, but not applicable for strong non-equilibrium states. A system may undergo a strong non-equilibrium initially and a weak non-equilibrium finally. It is quite difficult to find a universally valid way to decompose total concentrations, so constructing phase-concentration equations are more advisable to connect thermodynamic databases and phase-field simulations. In 2012, Steinbach et al. 103,104 proposed a phase-field model which directly gives the phase-concentration equations, but their equations cannot degenerate into Fick's diffusion equation from the interface to the bulk. More rigorous phaseconcentration equations are needed in the phase-field model for its universal applicability and convenient linkage with realistic thermodynamic databases.
The deposition reaction in charging process generally involves two phases (electrode and electrolyte). However, the intercalation reaction in discharging process may produce several new phases, i.e., three or more phases may be involved in the intercalation reaction. Presently, most phase-field models are rather successful in dealing with two-phase systems, but may be problematic when extended to multiphase systems (three phases and above). Nonrigorous mathematical treatment will be needed or nonphysical simulation results will appear; therefore, for phase-field simulations of RBs, further development is needed to build a highly rigorous and self-consistent phase-field model for multicomponent multiphase systems to describe, compare, analyze, predict, improve, and control properties of battery materials in the intercalation reaction. The expected phase-field model for multicomponent multiphase systems should have following properties: (1) all evolution equations are continuous in space and time, and the phaseconcentration equations should naturally degenerates into Fick's diffusion equations from the interface to the bulk; (2) the evolution equations of field variables directly reflect their constraint conditions and the same simulation result must be provided for different choice of independent variables; (3) the evolution equations of field variables can guarantee them in the proper range automatically, and truncation can be avoided in numerical simulations, except for the numerical errors; (4) field variables should be calculated strictly according to the evolution equation and free energy functional and non-logical mathematical skills must be removed from numerical simulations; and (5) it can be built into a sustainable open-source phase-field simulation community.
Problems to be considered. Intercalation of Li into the electrode leads to plastic deformation, which is related to the intercalation rate and viscoelastic behavior of the electrode, which is obviously different from the plastic deformation under external force. Subsequent phase-field simulations need to consider the quantitative relation between the above two factors and plastic deformation. The elasticity tensor is relevant to Li-ion concentration. At present, the function of elasticity tensor with respect to Li-ion concentration is assumed to be linear. However, this assumption is comparatively rough and a reasonable model is needed to reflect their actual relation in subsequent simulations involving stress. The plastic deformation resulting from cracks in the electrode will expose new surface of the electrode, where the electrochemical reaction and Li nucleation occur. Accordingly, the effect of plastic deformation on Li nucleation should be included in further simulations involving stress. The battery electrode is generally made of many grains mixed by adhesive. The shape of grains affects the Li-ion transport and stress distribution in the electrode 105,106 . This resulting stress sequentially affects the chemical potential, intercalation potential, deposition potential, and battery capacity. The grain-shape effect of the electrode should be considered in future simulations involving stress.
With the recent developments in solid-state electrolyte, the conductivity of the solid-state electrolyte at room temperature can be comparable or even larger than that of the liquid electrolyte. However, only few studies have been carried out to analyze the Li dendrite growth within solid electrolyte and liquid crystalline electrolyte. The nucleation and dendrite growth mechanism for solid and liquid crystalline electrolytes should be further studied. Based on the experimental studies about aqueous hybrid electrolyte battery in recent years, it has been found that hybrid ions can inhibit the dendrite growth 107,108 . However, the inhibition mechanism of dendrite growth in hybrid electrolyte is still unknown and needs deep investigations. Compared with ionic deposition in liquid electrolyte batteries, that in solid electrolyte batteries 109 depends on not only the current density, temperature, and electrolyte but also elasticity tensor of solid electrolyte, viscoplastic behavior of the electrode, and coating stress. These factors need to be considered in subsequent simulations 110 .

CONCLUDING REMARKS AND PERSPECTIVES
To our understanding, phase-field simulations can be classified into three levels. The first level describes microstructure evolution of materials, and compares microstructure evolution of similar materials to ascertain a material with better properties, which is a mixture of qualitative, semi-quantitative, and quantitative investigations. The second level analyzes the link between microstructure and performance, and predicts performance evolution of the system according to its simulation results, which is a mixture of semi-quantitative and quantitative investigations. The third level improves performance of the system based on its influence factors and controls its better service performance by adjusting its physical environment, which is a quantitative investigation. At present, phase-field simulations in RBs belong to the first two levels; specifically, most simulations are still in the first level. Progress of phase-field simulations for RBs toward quantification needs accurate thermodynamic and kinetic databases of the system. However, an RB is such a complicated system that many mechanical, chemical, thermal, and electric changes are involved in its charging or discharging processes. Therefore, the further phase-field model needs to combine DFT and MD to be developed into a multiscale model, whereby adjustable parameters can be eliminated and all parameters in the simulations can be obtained in terms of theoretical calculations or experimental measurements. After all, validity and rationality of the phase-field method, to a large extent, relies on physical, chemical, and mechanical parameters, and accurate thermodynamic and kinetic databases of the system.