Sodium Ion Batteries Particles: Phase-Field Modeling with Coupling of Cahn-Hilliard Equation and Finite Deformation Elasticity

The cathode material Na x FePO 4 of sodium-ion batteries shows phase changes during intercalation. In this work, a phase-ﬁeld model for Na x FePO 4 is studied for the ﬁrst time. The Cahn-Hilliard diffusion equation coupled to ﬁnite deformation elasticity is derived. Two ﬁnite deformation elasticity formulations based on elastic Green strain and logarithmic elastic strain, respectively, are compared. The material parameters for Na x FePO 4 are determined. We implemented the model in COMSOL Multiphysics for a spherically symmetric problem of sodium insertion into or extraction from a cathodic particle made of Na x FePO 4 . The model captures the important feature of phase segregation into a sodium-poor phase FePO 4 and a sodium-rich phase Na 2 / 3 FePO 4 . There is a visible difference for the concentration and stress between the small deformation theory and the ﬁnite deformation theories. Furthermore, we compare the two cathode materials Na x FePO 4 and Li x FePO 4 of lithium-ion batteries to each other in terms of phase changes and stresses, and show that although the miscibility gap of Na x FePO 4 is smaller than that of Li x FePO 4 , the stresses in the cathode material Na x FePO 4 are higher in the phase segregated state. As a result, the suppression of phase segregation by the elastic strain energy is more easily achieved in Na x FePO 4 compared to Li x FePO 4 .

Lithium-ion batteries (LIBs) have attracted intensive attention for electrochemical energy storage over the past decades. The massive use of LIBs combined with the limited and unevenly distributed lithium source will drive up the prices of lithium, and high cost remains a critical problem for the development of LIBs. 1 In contrast to lithium, the wide availability, abundance and low cost of sodium on earth make sodium-ion batteries (NIBs) suitable for large scale energy storage devices in which high energy density becomes less critical. 2 For comparison, properties of lithium and sodium are summarized in Table  I. Recently, NIBs have been considered as a promising alternative to LIBs, since sodium and lithium exhibit similar chemical properties so that sodium chemistry could be applied to a similar battery system, and the fundamental principles of the NIBs and LIBs are identical. 2 Just like LIBs, the electrochemical processes in the electrodes of NIBs are also coupled to mechanical properties. Insertion of sodium into the host material leads to strains, and as a consequence, stresses are induced. The stresses associated with the related volume changes may lead to fracture of electrode particles. For example, in LIB systems a volume expansion of more than 300% of silicon particles has been observed, 3 which can cause particle fracture and capacity loss. On the other hand, for thermodynamical reasons, there is a contribution of the stresses to the driving force for diffusion in the host material.
Among all cathode materials for NIBs, phosphate based cathode materials are among the best candidates, due to their thermal stability and higher voltage. 1 Olivine N a x FeP O 4 (NaFPO) has the highest theoretical specific capacity (154m Ahg −1 ) compared to the other phosphate polyanion cathode materials (N aV P O 4 F, N a 3 V 2 (P O 4 ) 2 F 3 and N a 2 FeP O 4 F, etc.), 1 which makes it a promising cathode material for NIBs. Similar to olivine Li x FeP O 4 (LiFPO), olivine NaFPO also shows phase changes during sodium insertion or extraction. 2,[4][5][6][7] For the modeling of phase changes, Huttin 8 has formulated a phasefield model employing the Cahn-Hilliard equation considering the cathode material Li x Mn 2 O 4 (LMO) of LIBs. Phase-field simulations based on the Cahn-Hilliard equation for LiFPO particles of LIBs are also studied in Refs. 9,10. A thermodynamic phase-field model based on the Cahn-Hilliard equation relies on a continuous order parameter which is a conserved quantity, 11 thus, leading to diffuse interfaces between adjacent phases with no need for the cumbersome tracking of the position of a sharp interface. Our recent work 12 has discussed z E-mail: tao.zhang@kit.edu the nonlocality of the Cahn-Hilliard theory, and demonstrated that it is weakly nonlocal. In this work, a phase-field model based on the Cahn-Hilliard equation for NaFPO of NIBs will be studied for the first time, and the relevant material parameters of NaFPO for the phasefield model will be determined. As far as we know, no work has been published for the phase-field modeling of the cathode materials of NIBs by now.
As for the coupling between diffusion and mechanics, Huttin and Kamlah 13 considered the coupling between the Cahn-Hilliard equation and small deformation theory (SDT) for spherical particles of LMO. A general phase-field theory for the coupling of the Cahn-Hilliard equation and finite deformation elasticity based on the local dissipation inequality and the so-called local microforce balance has been presented by Gurtin. 14 In this theory, the external microforce itself remains a quantity free to choose, and by taking it equal to zero, the classical Cahn-Hilliard equation is obtained. Walk et al. 15 have identified the influence of the external microforce in the chemical potential for spherical LMO particles.
Anand 16 derived a theory for lithium insertion modeled by the Cahn-Hilliard equation combined with finite deformation elastoplasticity based on the multiplicative decomposition of the deformation gradient, in which the external microforce is identified as related to the Mandel stress in the chemical potential. To this end, the socalled macroforce and microforce balances are evaluated by following the virtual-power approach of Germain 17 and Gurtin. 18 Subsequently, Di Leo et al. 19 formulated a continuum model which coupled the Cahn-Hilliard-type phase-field theory with finite deformation elasticity based on logarithmic elastic strain. However, Anand's model 16,19 ignored the fact that logarithmic elastic strain is dependent on the species concentration in the stress chemical potential, although logarithmic elastic strain is also a constitutive variable. As a result, the extremely complicated and important term induced by the concentration dependence of the logarithmic elastic strain in the stress chemical potential is not accounted for. The concentration dependence of the strain tensor in the chemical potential is also confirmed by Gurtin. 14 Interestingly, the external microforce term related to the Mandel stress in Anand's model 16,19 is just equivalent to the ignored term induced by the concentration dependence of the logarithmic elastic strain, which also is verified by Refs. 15 and 20. Walk et al. 15 reported that for SDT the identification of the external microforce in Anand's model 16,19 leads to a doubling of the mechanical coupling term in the chemical potential. Zhao et al. 20 derived a Cahn-Hilliard phase-field model coupled to mechanics based on the Neo-Hookean elasticity without introduction of the external microforce in the chemical potential, and show that their model agrees with Anand's model 16,19 which ignored the concentration dependence of the strain tensor in the chemical potential. As remarked by Gurtin,14 no constitutive relation is specified for the external force. Rather, specifying the external microforce by invoking the so-called principle of virtual power may be considered an option. However, in view of the above discussion we do not find it appropriate to adopt this here. In our paper, we will formulate a coupled phase-field theory for NaFPO modeled by the Cahn-Hilliard theory coupled to finite deformation elasticity formulations, and the concentration dependence of the strain tensor in the chemical potential will be considered. For the mechanical part, two different finite deformation elasticity formulations will be introduced and compared.
The paper is organized as follows. In Theory section, we derive the coupled model of the Cahn-Hilliard diffusion theory and finite deformation elasticity, including two mentioned mechanical theories based on elastic Green strain and logarithmic elastic strain, respectively. The governing equations and boundary conditions are summarized. The resulting set of equations are implemented in COMSOL Multiphysics for the solution of the initial-boundary-value problem by the finite element method. We describe the determination of the material parameters for NaFPO cathode material in Material Parameters section. We present and discuss phase changes and stresses in spherical NaFPO electrode particles, and compare the cathode material NaFPO of NIBs to LiFPO of LIBs in terms of phase changes and stresses in Results and discussion section. Finally, we conclude our study in Conclusions section.

Theory
In this section, we will first motivate the phase-field model for the pure diffusion problem. Then, the coupling to mechanics will be summarized. For the mechanical part, we will introduce two different finite deformation theories including elasticity based on elastic Green strain and elasticity based on logarithmic elastic strain.

Cahn-Hilliard model for diffusion.-
The description of diffusion and phase changes within an electrode material can be modeled by the Cahn-Hilliard theory. 11 In the phase-field model, we introduce as an order parameter, depending continuously on space, the species concentration c, which is measured in mol per unit volume. The free energy density consisting of two parts is given by ψ(c, gradc) = ψ mwp (c) + ψ gd (gradc), [1] wherec is the dimensionless concentration, which is normalized with respect to the maximum species concentration c max asc = c/c max . The homogeneous free energy density ψ mwp is a multiwell potential defining the respective phases. Based on the References 13,24 the homogeneous free energy density exhibits a double-well structure. The existence of this zone of concavity indicates that homogeneous species concentration states do not always ensure the system free energy to be minimal. In the concavity zone, the system becomes unstable toward phase segregation. The Maxwell construction, which connects the neighborhoods of the two minima by a common tangent, predicts that the system splits into a two phase system and the chemical potentials of both the species-poor phase and species-rich phase are equal.
In this sense, the homogeneous free energy density is a multiwell potential where the wells define the respective phases. The second term on the right hand side of Equation 1 represents the gradient energy leading to a diffuse interface between two adjacent phases, which is given by Here k B is the Boltzmann constant, N A is the Avogadro constant and T re f is the reference temperature. Furthermore, λ is a material constant with units of length squared controlling the thickness of the diffuse interface and | · | denotes the norm of a vector.
The system free energy of some domain of volume V is According to Cahn and Hilliard, 11 the system free energy (c) can be interpreted physically such that, to the first approximation, the free energy of a small volume of nonuniform solution can be expressed as the sum of two terms, one being the free energy that this volume would have in a homogeneous solution and the other a "gradient energy". The homogeneous free energy density for a two phase material can be expressed as 13,24 ψ mwp (c) ) . [4] The first two terms on the right hand side of Equation 4 represent the interaction energy, where positive values of α 1 characterize the energy of inserting a species into the host material, and negative values of α 2 indicate the interaction of neighboring species to be attractive. Concerning α 1 , it is hard to relate a physical meaning to non-positive values of α 1 . On the other hand, it has to be noted that α 1 does not occur in the Cahn-Hilliard evolution equation for concentration. According to Huttin,8 there is a critical temperature, namely T c = −1/4α 2 T re f , above which the homogeneous free energy density is of convex shape, representing an ideal solution. For temperatures T below the critical temperature, the free energy density exhibits a zone of concavity where the homogeneous concentration states are not stable states of the system, and phase segregation becomes possible. The zone of concavity corresponds to the condition that the second order derivative of the homogeneous free energy density with respect to the concentration is negative. This inequality is never fulfilled if α 2 is not negative. Thus, for a system of noninteracting species, i.e. α 2 = 0, or for a system of species that repel each other, i.e. α 2 > 0, the homogeneous free energy density keeps a convex shape for all T > 0. For an attractively interacting species system, i.e. α 2 < 0, depending on the temperature, the homogeneous free energy density may exhibit a zone of concavity. Therefore, at T = T re f , the attraction is strong enough to initiate phase segregation for α 2 < −4. The terms multiplied by absolute temperature T represent the entropy of mixing, wherec lnc represents the entropy of mixing valid at low concentration, and (1 −c) ln (1 −c) describes the non-dilute solution behavior and represents the entropy of mixing responsible for the saturation effect. 20,25 Journal of The Electrochemical Society, 165 (10) A1997-A2007 (2018) A1999 Table II. Basic fields in the kinematics.

Descriptions
Basic fields The driving force for diffusion is expressed as the gradient of the chemical potential. Based on the local dissipation inequality and the so-called local microforce balance, 14 the chemical potential can be obtained as The external microforce γ is a contribution in the so-called local microforce balance, which can be chosen arbitrarily. 14,15 In our case, the external microforce is neglected (γ = 0). Gurtin 14 also obtains the constitutive equation for the mass flux as the Onsager relation where the mobility tensor M is non-negative definite. The diffusion equation governing species transport is based on the conservation of mass, namelyċ = −div( J ). [7] Combining Equations 6 and 7 yields the traditional Cahn-Hilliard equationċ = div (M(c)grad μ) , [8] where in this special case the isotropic mobility M(c) is a non-negative function. Table II.

Finite deformation elasticity.-Kinematics.-The kinematics developed here is related to the basic fields shown in
Consider a motion x = χ( X , t) of a homogeneous isotropic body B, where X denotes an arbitrary material point of B in the reference configuration, and x is the spatial position of this material point in the current configuration. The deformation gradient is then defined as F = Grad χ, [9] where Grad is the gradient operator calculated with respect to the reference configuration. The kinematical basis for the coupling between finite deformation elasticity and diffusion is the multiplicative decomposition 16 as illustrated in Fig. 1. Here F e is the elastic deformation gradient that is a mapping from the intermediate configuration onto the current configuration, and the concentration induced deformation gradient F s caused by species insertion or extraction is a mapping from the reference configuration onto the intermediate configuration. F s is assumed to be purely volumetric in the form [11] where c R is the species concentration, measured in mol per unit reference volume, is the partial molar volume, and I is the second order unit tensor. The volume ratio J between a current volume element dV c and an initial volume element dV i during species insertion or extraction is defined as Using Equation 10 yields where J e = detF e and J s = detF s . The right and left polar decompositions of F e are expressed by where R e is a rotation, and U e and V e are the symmetric, positivedefinite elastic right and left stretch tensors: In addition, the right elastic Cauchy-Green tensor is The elastic Green strain tensor is then defined as Next, the spectral representations of U e and V e are , [18] respectively where ( r 1 , r 2 , r 3 ) and (λ e 1 , λ e 2 , λ e 3 ) are the orthonormal eigenvectors and the positive eigenvalues of U e respectively. According to Anand, 16 the logarithmic elastic strain tensor and the spatial logarithmic elastic strain tensor can be obtained as Journal of The Electrochemical Society, 165 (10) A1997-A2007 (2018) Elasticity based on elastic Green strain.-According to Gurtin, 14 the thermodynamical potential relation for the first Piola-Kirchhoff stress T R is derived from the local dissipation inequality as Here, the first Piola-Kirchhoff stress T R is related to the symmetric Cauchy stress T in the current configuration by The second Piola-Kirchhoff stressT e related to intermediate configuration is symmetric considering the symmetry of the Cauchy stress [23] Using Equations 21-23, we can obtain the relatioñ The free energy density of the coupled theory, now per unit reference volume, includes three parts: Here ψ mwp R (c) and ψ gd R have been previously discussed as the multiwell potential and the gradient energy. The additional term ψ cp R (E e , c R ) is the coupling energy in terms of elastic Green strain defining the coupling between diffusion and mechanics, which is assumed to be given by Here C is the elasticity tensor which is taken to be constant and isotropic. Accordingly, G is the shear modulus, and ν is the Poisson number. Using Equations 25 and 26, Equation 24 then yields the isotropic elasticity law based on elastic Green strain asT [28] Elasticity based on logarithmic elastic strain.-In the following, we will introduce another isotropic elasticity law based on logarithmic elastic strain. According to Anand, 16 the second Piola-Kirchhoff stress T e is derived from the free energy density according tõ Here [33] The free energy density of the coupled theory in this case is while the coupling energy in terms of logarithmic elastic strain reads as [35] Using Equations 34 and 35, Equation 32 then yields the isotropic elasticity law based on logarithmic elastic strain as [36] The coupled theory.-Using Equation 5, the chemical potential considering the coupling between diffusion and mechanics is given by where μ cp R is the coupling chemical potential. We consider two forms of it for the two different finite deformation elasticities: [39] Here, the term multiplied by ∂E e ∂c R or ∂E e log ∂c R represents the concentration dependence of the strain tensor in the chemical potential. It should be noticed that this term is equivalent to the external microforce term related to the Mandel stress in Anand's model. 16,19 We choose an isotropic mobility according to with the function which is symmetric in the range between zero and maximum concentration and in which D 0 is the diffusion coefficient. Finally, based on the balances of mass and momentum, respectively, we can obtain the field equations Div T R = 0. [43] Combined with the constitutive equations introduced above, the field equations can be taken as partial differential equations for concentration c R and displacement vector u, which need to be solved for given initial and boundary conditions. This is a fourth-order nonlinear initial-boundary-value problem.
In small deformation theory, the linear strain tensor is ε i j = 1/2 u i, j + u j,i with u i being the displacement vector, and Hooke's law gives for the stress relation T = C : (ε − ε s ), where ε s = 1/3 (c −c 0 )I is the stress-free strain induced by species insertion or extraction. Boundary and initial conditions.-A spherical particle of radius R 0 under spherically symmetric boundary conditions is considered. Here, spherical coordinates are introduced and all fields are expressed as a function of the time t and the radial coordinate 0 ≤ R ≤ R 0 of the reference configuration in the form [44] [45] The boundary conditions are sketched in Fig. 2. At the surface, the particle is assumed to be stress free, where n R refers to the outgoing unit vector normal to the particle surface. We choose a spatially independent mass flux at the surface as [47] Here C is the C-rate, and C = n means that the amount of a species of a fully charged particle would flow into the particle within 1/n hours. Once the maximum concentration c max is reached at the surface, the mass flux will be stopped. We also impose the "natural" or "variational" boundary condition at the surface for the fourth order Cahn-Hilliard equation. 10 Here neglecting surface wetting, the "natural" boundary condition is expressed as When the interface between phases meets the particle surface, Equation 48 enforces that it is perpendicular to the surface. 20 In addition, at the particle center, the boundary conditions are to be satisfied. Here, Equation 49 guarantees the continuity of the particle at the center, and Equations 49-51 are needed to ensure the spherical symmetry. What is more, these three boundary conditions ensure the physical requirement that the flux at the particle center vanishes. 8 Finally, the initial conditions are given by c R (R, 0) = 0. [53]

Material Parameters
As mentioned before, we simulate particles of the cathode material sodium iron phosphate NaFPO. We consider a typical active nanoparticle of size R 0 = 150 nm for the particle radius. According to Nakayama et al., 26 the diffusion coefficient of Na-ions in NaFPO is an order of magnitude lower than that of Li-ions in LiFPO, which also can be confirmed by comparison to Zhu et al. 4 The diffusion coefficient of LiFPO is 1 × 10 −14 m 2 /s, 10,27 therefore, we use the value D 0 = 1 × 10 −15 m 2 /s as the diffusion coefficient of NaFPO. In the following, we will specify the values of the other parameters for NaFPO including c max , α 1 , α 2 ,λ, partial molar volume and Young's modulus E 0 .
Determination of c max .-The maximum concentration c max with units of mol/m 3 means the maximum content of a certain species the host material can accept, which can be expressed by 8 where V 0 is the volume occupied by one species atom. This can be obtained by with n being the number of atoms of the species that the host material can accept in a unit cell of volume V cell . According to Ref. 28, the structures of olivine LiFPO and olivine NaFPO are the same. Based on the crystal structure of LiFPO as shown in Ref. 29 and the lattice parameters in Table III, we can calculate the maximum lithium concentration of LiFPO by Equation 54, which gives 2.29 × 10 4 mol/m 3 , matching the value used by Zeng and Bazant, 10 and Delacourt and Safari. 30 Therefore, in the same way, the maximum sodium concentration of NaFPO can be obtained as c max = 2.1 × 10 4 mol/m 3 , which is smaller than that of LiFPO due to the larger unit cell volume for NaFPO.
Determination of α 1 and α 2 .-According to Lu et al., 5 at room temperature, the phase diagram of olivine NaFPO consists of two regions. For 0 < x < 2/3, phase segregation into a sodium-poor phase FeP O 4 and a sodium-rich phase N a 2/3 FeP O 4 is found to be favorable, which means that this is a two-phase region. (Following common praxis, we replace in chemical formulas dimensionless concentration c by x.) Based on the unit cell parameters and volume as functions of x in NaFPO as shown in Ref. 5, two different values for the unit cell parameters and volume, respectively, can be identified when x increases to around 0.08, which is the manifestation of the occurrence of the phase segregation. Therefore, we can obtain the normalized sodium concentration for the initiation of phase segregation, which is around 0.08. In contrast to LiFPO  With regard to the homogeneous free energy density in the twophase region of NaFPO, we assume This is the classical homogeneous free energy density function for a two phase material, 8,24 which we have formulated such that it is limited to the range 0 < x < 2/3. Consequently, we fit the homogeneous free energy density (56) to the above experimental data 5 with respect to the parameters α 1 and α 2 , as shown in Fig. 3. The best fit was achieved with α 1 = 5 and α 2 = −15 at room temperature. Fig. 3 shows the fitting result with the normalized homogeneous free energy density in the two-phase regionψ mwp plotted versus normalized concentrationc. It can be seen thatψ mwp exhibits a doublewell structure, such that two different relative minima A and B occur, characterizing two phases of different equilibrium values for the sodium concentration. The two minima A and B correspond to the two phases FeP O 4 and N a 2/3 FeP O 4 , respectively. The two zones AC and BD are the "nucleation zones", and phase segregation occurs upon sufficient disturbance of the system's equilibrium. In the inner zone of concavity between the two points of inflection C and D, which is denoted the "spinodal decomposition zone", homogeneous sodium concentration states are unstable and phase segregation is initiated in any case. The Maxwell construction represents the volume fractions of a phase segregated system. The inflection points correspond to the second order derivative of ψ mwp with repect toc, as shown in Fig. 4, and the inflection point C matches the experimental data from Ref.

5.
Determination of λ.-The so-called reference value λ 0 for λ can be estimated as 8 and the relation between λ 0 and λ is where α is constant. Using Equation 57, we obtain λ 0 = 9.1 × 10 −19 m 2 for NaFPO. In absence of any further information, we pick α = 20 as the value for NaFPO which is the one digit number in the order of magnitude of the values for LiFPO and LMO, see Table IV. Therefore, we get   Table IV, as well.
Determination of .-The partial molar volume plays a role analogous to a thermal expansion coefficient, meaning we can calculate by the relation ε s = 1/3 (c −c 0 )I. Thus, based on Table V, the partial molar volume of LiFPO is obtained as which matches the value used by Song et al. 25 Therefore, in the same way, the partial molar volume of NaFPO can be calculated as This value of is larger than that of LiFPO, which is consistent with the fact that sodium has a larger cation radius than lithium as shown in Table I. It should be noticed that we do not consider the possibility of a concentration dependence of .  32 Based on the above data, due to the larger ionic radius of Na-ion compared to that of Li-ion, we regard the two digit estimate 120 G Pa as Young's modulus of NaFPO.
The material parameters for the two cathode materials LiFPO and NaFPO are summarized in Table VI. The resulting set of equations has been implemented in COMSOL Multiphysics for solution by the finite element method.

Results and Discussion
In this section, we will consider the quasistatic insertion and extraction of sodium for a particle of NaFPO at C = 0.001. In our simulations, we exclusively focus on the two-phase region of NaFPO (0 < x < 2/3). In the figures, c, c R and T H , T H R denote concentration and hydrostatic stress of small and finite deformation elasticity, respectively, where T H = 1/3T ii is the hydrostatic stress in terms of the Cauchy stress. We define the normalized quantitȳ The average concentration c avg , also called "state of charge" (SOC) is Fig. 5, the markers represent values of the normalized average system free energȳ

Cahn-Hilliard model without mechanics.-In
during insertion and extraction, respectively, as function of c avg . For demonstration purposes, the plot of the homogeneous free energy density vs. dimensionless concentration is entered, as well. Markers on the homogeneous free energy density curve correspond to homogeneous states whereas markers nearby the path of the Maxwell construction correspond to phase segregated states, as they are illustrated in Figs. 6 and 7.    6 shows the normalized sodium concentration along the radial coordinate at different time instants during sodium insertion. We find that at the beginning of sodium insertion the concentration is homogeneous. Once c avg gets close to 0.08, which corresponds to the inflection point C in Fig. 3, phase segregation is initiated. Two different phases can be recognized in the particle, namely the sodium-poor phase FeP O 4 in the center corresponding to the first minimum A and the sodium-rich phase N a 2/3 FeP O 4 at the outside corresponding to the second minimum B. A smooth but very narrow interface with concentration changing rapidly but continuously separates them. When c avg grows up to 2/3, the intermediate phase N a 2/3 FeP O 4 occupies all of the particle.
Next, we consider the sodium extraction case shown in Fig. 7. At the beginning the system is in a homogeneous state and the corresponding concentration is 2/3. Once c avg is reduced to around 0.6, which corresponds to the inflection point D in Fig. 3, phase segregation is initiated during sodium extraction. In contrast to the insertion case, the sodium-poor phase FeP O 4 is at the outside while the sodium-rich phase N a 2/3 FeP O 4 is in the center. At the end of sodium extraction, the sodium-rich phase N a 2/3 FeP O 4 vanishes. It should be noticed that although to the particle surface a mass flux is applied, the concentration on the particle surface nearly stays at a constant value. This is due to the fact that because of the extremely low C-rate of C = 0.001  for quasistatic insertion and extraction of sodium, the system moves along a path of relaxed quasiequilibrium states.

Cahn-Hilliard model with mechanics.-Concerning mechanics,
we first study the coupling between the Cahn-Hilliard theory and small deformation theory (SDT), and then compare the different small and finite deformation theories. In order to study the coupling effect through the coupling energy on sodium diffusion and stress in the particle, we introduce a normalized Young's modulus as where E 0 is the value of Young's modulus shown in Table VI. Hilliard model with SDT.-Figs. 8 and 9 show normalized sodium concentration and hydrostatic stress along the radial coordinate at a stage of average concentration of c avg = 0.5 during sodium insertion based on SDT for differentĒ. We find that the difference in concentration between the two phases is reduced asĒ increases, which means that the solid solution limits of FeP O 4 and N a x FeP O 4 are gradually extended into the range of phase segregated states. In other words the coupling effect reduces the miscibility gap. Correspondingly, the hydrostatic stress magnitudes increase asĒ increases. The hydrostatic stress is constant and tensile in the low concentration  phase, and it drops rapidly across the interface to become compressive and constant in the high concentration phase. This change of sign is a result of mechanical equilibrium between the inner core and the outer shell of the particle. Note that on the other side forĒ = 1, the system is homogeneous and stress-free. Even though the average concentration value lies in the spinodal decomposition zone ofψ mwp , phase segregation does not occur. As illustrated in Fig. 10, for a high value of normalized Young's modulus, the cost of the coupling system energy¯ cp at a phase segregated state would be too high so that as a result of the competition between the different contributions to the system free energy the system stays in the homogeneous state in order to minimize the total system free energy. Therefore, the presence of a coupling energy in the total system free energy can lead to suppression of phase segregation.   Fig. 11. It can be recognized that there is a visible difference for the sodium concentration plots of SDT and the finite deformation elasticity formulations, with the difference between the two phases being less pronounced for SDT, while the sodium concentration plots of two finite deformation theories are almost the same. Compared to the experimental result, 7 we can find that the sodium concentration in the high concentration phase obtained from two finite deformation theories matches the experimental value for N a 0.6 FeP O 4 but there is a clear deviation from the result by SDT. On the other hand, the sodium concentration in the low concentration phase obtained from all three mechanics theories is always slightly below the experimental value for N a 0.08 FeP O 4 . In Fig. 12, the hydrostatic stress plots are clearly different for SDT on the one side and the two finite deformation theories on the other. There are higher stresses in the two phases for SDT. This is attributed to the fact that SDT does not account for the volume swelling of the particle during sodium insertion. Comparing the two Figure 13. Normalized sodium concentration versus normalized radial coordinate for variousĒ during sodium insertion based on different mechanics theories. finite deformation theories to each other, it is found that there is a negligible difference in stress between them.

Comparison of different mechanics
As mentioned before, the contribution of the coupling energy to the total system free energy can lead to suppression of phase segregation. Here, we discuss the critical value of normalized Young's modulusĒ c of NaFPO above which phase segregation cannot arise in the particle. Fig. 13 shows the normalized sodium concentration along the radial coordinate for various values ofĒ during sodium insertion based on different mechanics theories. For SDT, it can be seen that phase segregation can arise at a stage of average concentration of c avg = 0.33 whenĒ = 0.385 but sodium concentration is homogeneous over the particle whenĒ = 0.386, soĒ c is 0.385 for SDT. In the same way, we findĒ c = 0.409 for the two finite deformation elasticity formulations, which is a little larger than that for SDT.
Comparison of different cathode materials.-We now compare the two cathode materials NaFPO and LiFPO in terms of phase changes and hydrostatic stress. Figs. 14 and 15 show normalized concentration and hydrostatic stress along the radial coordinate for various particle radii R 0 at a stage of average concentration of c avg = 0.5 for the two cathode materials during insertion based on different finite deformation theories. In Fig. 14, we recognize that there is an obviously reduced miscibility gap for a NaFPO particle compared to  a LiFPO particle. This is consistent with the experimental fact that NaFPO has a two-phase region in the range 0 < x < 2/3 while phase segregation for LiFPO occurs in the range 0 < x < 1. 5 Additionally, the interface location in a NaFPO particle is more close to the particle center, and, furthermore, the concentration in the two phases is the same for different particle radii.
Although the miscibility gap of NaFPO is much smaller than that of LiFPO, the stresses in a NaFPO particle shown in Fig. 15 are larger compared to those in a LiFPO particle. This can be due to a larger expansion during phase changes for NaFPO as shown in Fig.  16 describing the plots of volume ratio J at the particle surface. Here, it should be mentioned again that our simulations only consider the two-phase region of NaFPO (0 < x < 2/3). We can find that the volume expansion from FeP O 4 to N a 2/3 FeP O 4 is quite large, namely about 12.8%, which is nearly 2 times that for LiFPO changing from FeP O 4 to Li Fe P O 4 . This is consistent with the reports. [5][6][7] Correspondingly, as shown in Fig. 17, the maximum hydrostatic stress magnitude in a NaFPO particle during phase changes is always larger than that in a LiFPO particle. The larger stresses in a NaFPO particle may explain the existence of a wide range of solid solution N a x FeP O 4 (2/3 < x < 1) to avoid even higher stresses in the NaFPO particle. Casas-Cabanas et al. 6 conclude that the larger stresses in a NaFPO particle can be the explanation for the existence of an intermediate phase. Indeed, the formation of an intermediate phase acts as a buffer between FeP O 4 and N a FeP O 4 providing elasticity to the structure. On the other hand, the critical value of normalized Young's modulusĒ c of LiFPO for surpression of phase segregation is equal to 1 (we do not show this result here), which is larger than that of NaFPO, as shown before. Therefore, it is easier for NaFPO to reach the massive coupling energy to suppress phase segregation compared to LiFPO due to the existence of larger stresses in a NaFPO particle.
As can be seen in Fig. 17, the LiFPO particle surface expands more seriously around c avg = 0.15. This is due to phase segregation being already initiated for LiFPO but not yet for NaFPO which is still in the homogeneous state. As a result, the maximum hydrostatic stress magnitude in a LiFPO particle is far greater than that in a NaFPO particle at this insertion state. At c avg = 2/3, the maximum hydrostatic stress magnitude in a NaFPO particle is close to zero while the LiFPO particle still displays high stresses. This is attributed to the intermediate phase N a 2/3 FeP O 4 occupying then the whole particle while LiFPO is still in a phase segregated state.

Conclusions
A phase-field model for the cathode material NaFPO of NIBs is studied for the first time. A coupled phase-field model for Cahn-Hilliard theory and elasticity has been derived. For the mechanical part, besides small deformation theory (SDT) two different finite deformation elasticity formulations are introduced and compared. We have taken into account the concentration dependence of the strain tensor in the chemical potential which seems to have been ignored by other researchers. As a major novelty, the material parameters for NaFPO are determined. For example, α 1 , α 2 and λ, all of which are the key parameters in the phase-field model, are determined. The determination of these key parameters provides a significant input for the future phase-field work for NaFPO. We implemented the fourth-order nonlinear initial-boundary-value problem of the model in COMSOL Multiphysics to solve by the finite element method the spherically symmetric problem of sodium insertion into or extraction from a NaFPO particle of NIBs.
Our model captures the important feature that distinguishes NaFPO from LiFPO, i.e, phase segregation into a sodium-poor phase FeP O 4 and a sodium-rich phase N a 2/3 FeP O 4 . The difference for the concentration and stress plots between SDT and the finite deformation theories cannot be neglected for NaFPO. In particular, the stresses in the two phases are higher for SDT. The above difference suggests that the finite deformation elasticity is preferred for the cathode material NaFPO. This is different from that the small deformation theory has a sufficient capacity to represent the deformation of LiFPO. Although the results from the two finite deformation theories are almost the same, according to our experience, elasticity based on logarithmic elastic strain is numerically more efficient than elasticity based on elastic Green strain. On the other hand, we find that the miscibility gap of NaFPO is smaller than that of LiFPO, but the stresses in a NaFPO particle during phase segregation are larger, which may explain the existence of a wide range of solid solution N a x FeP O 4 for 2/3 < x < 1 to avoid even higher stresses in the NaFPO particle. In addition, the suppression of phase segregation by the elastic strain energy is more pronounced in NaFPO compared to LiFPO.
In this work, we just focus on phase changes in the two-phase region of NaFPO (0 < x < 2/3), ignoring the single-phase region (2/3 < x < 1). We will extend the current model into the whole region (0 < x < 1) in our future work. As both the partial molar volume and Young's modulus are regarded as constants here, another major subject of future study is the effect of the concentration dependence of these two quantities for NaFPO.