Quantitative Phase Field Model for Electrochemical Systems

Modeling microstructure evolution in electrochemical systems is vital for understanding the mechanism of various electrochemical processes. In this work, we propose a general phase field framework that is fully variational and thus guarantees that the energy decreases upon evolution in an isothermal system. The bulk and interface free energies are decoupled using a grand potential formulation to enhance numerical efficiency. The variational definition of the overpotential is used, and the reaction kinetics is incorporated into the evolution equation for the phase field to correctly capture capillary effects and eliminate additional model parameter calibrations. A higher-order kinetic correction is derived to accurately reproduce general reaction models such as the Butler-Volmer, Marcus, and Marcus-Hush-Chidsey models. Electrostatic potentials in the electrode and the electrolyte are considered separately as independent variables, providing additional freedom to capture the interfacial potential jump. To handle realistic materials and processing parameters for practical applications, a driving force extension method is used to enhance the grid size by three orders of magnitude. Finally, we comprehensively verify our phase field model using classical electrochemical theory.


Introduction
Electrochemical systems have many practical applications, such as electrodeposition and electrochemical energy storage.It is believed that electrochemical systems, like batteries, play an essential role in the global effort to achieve carbon neutrality.Understanding detailed electrochemical processes is key for performance optimization, materials and configurations design, and charging and discharging optimization.However, electrochemical systems are heterogeneous and composed of electrodes and electrolytes.Microstructural evolution at the mesoscale, like dendrite growth and isolated-metallic particle formation in batteries, is common in electrochemical processes.Thus, tracking the microstructure evolution is essential for modeling electrochemical processes in electrochemical systems at the mesoscale.
The phase field method has become an increasingly popular tool for interface tracking due to its flexibility in handling complex morphologies and coupled physical processes [1].However, there has been a historical challenge to develop practical and useful phase field models of electrochemical processes, which has led to a wide variety of proposed models.
Guyer et al. [2,3] developed a variational and thermodynamically consistent model for electrodeposition with linear kinetics that resolves the detailed structure of the electrical double layer.Angstrom scale grid size is needed for this model, so simulations beyond one-dimensional domains are rare.Shibuta et al. [4,5] and Okajima et al. [6] developed a variational model based on the classical Kim-Kim-Suzuki model of solidification [7] but did not explicitly model the electrical double layer.This model allows for a much larger diffuse interface width due to decoupling the bulk and the surface energies.However, a local equilibrium constraint must be fulfilled at each grid point, leading to a high computational cost for general free energy formulations.Cogswell [8] used the grand potential energy to eliminate explicit handling of the local equilibrium constraint.Besides increasing the diffuse interface width, one key advance in the thermodynamically consistent phase field modeling is that the overpotential arises from the variational principle, first given by Bazant [9].Many models employ this definition of the overpotential [10,8].Moreover, the general free energy form given by García et.al. [11] and related variational principles is another important step for phase field models of electrochemical systems.
Another important aspect in quantitative phase field modeling of electrochemical processes is the need to incorporate experimentally accessible reaction kinetics.There are generally two ways to include the reaction kinetics in the governing equations.The first introduces the reaction kinetics as a source term, like in the Allen-Cahn and the Cahn-Hilliard reaction models [9,12,13,14,15,16,17].These models enable simulations of a variety of phenomena in electrochemical systems, e.g., as demonstrated by García's group [14,17].However, this approach is nonvariational and thus does not guarantee that the energy will decrease during evolution.Therefore, the resulting model must be carefully verified, especially for the reaction kinetics and the capillary overpotential.Specifically, additional model parameter (phase field mobility) calibration is needed that may not be required in the standard phase field models in order to capture the correct reaction kinetics [13].A second way is to include the reaction kinetics in the evolution equation for the phase field variable [8,10].The source terms come out naturally from the variational derivation of the governing equations in a thermodynamically consistent way.Note that in the work by Chen et.al. [10] and many related works, e.g., [18,19,20], the phase field equation is reformulated into an Allen-Cahn reaction equation, so it is not clear if realistic capillary effects are captured and model parameter calibration (L η in [10]) is still needed.In addition, direct incorporation of the reaction kinetics into the phase field equation may lead to deviation of the phase field kinetics from the desired kinetics at a large overpotential.
Another aspect that requires careful treatment in developing a phase field model is the potential jump at the electrode-electrolyte interface.The electrostatic potential typically exhibits a large jump (can be more than 1 V) across the interface due to the existence of the electrical double layer.However, except for models that explicitly resolve the electrical double layer [2,3], the potential jump is frequently overlooked, or a constant jump is assumed, and a single electrostatic potential that continuously varies across the interface is used [12,10,8,18,14].
Finally, to guarantee that the phase field model can provide quantitative predictions, it is important to test whether it can reproduce the existing results of classical electrochemical models.In addition, these tests should be performed with experimentally accessible parameters on realistic length and time scales.However, a comprehensive model verification was not found in the literature.
In this work, we propose a general phase field framework that is fully variational (guarantees energy decrease) and is thermodynamically consistent.We propose a grand potential formulation to decouple the bulk and the surface free energies and eliminate explicitly solving conditions.E cell = −V app is the half-cell voltage.Illustration of field variables to describe the status of the system: (c) the phase field variable ξ varying across a length scale called the diffuse interface width l, (d) the electrostatic potential ϕ and the virtual electrostatic potentials ϕ s and ϕ l for each phase, and (e) the concentration c i of species i and the virtual concentrations c s i and c l i for each phase.
the local equilibrium constraint at each grid point.The classical local equilibrium constraint [21] is extended using the diffusional electrochemical potential.The variational definition of the overpotential follows from our approach, and the reaction kinetics is considered in the evolution equation for the phase field variable to avoid manual model parameter calibration.A higher-order kinetic correction is derived to reproduce general kinetic models for any overpotential accurately.We consider both electrostatic potentials in the electrode and the electrolyte as independent variables to provide additional freedom to capture the potential jump across the electrode-electrolyte interface.Potential jump varying along the interface can be correctly captured.To handle realistic materials and processing parameters with a micrometer-scale grid size, we apply a driving force extension method [22] to increase the grid size by three orders of magnitude for practical applications.Finally, we provide a comprehensive verification of our model, including the equilibrium behavior, the role of interfacial energy or capillarity, the ability to recover general reaction kinetics, the coupling between diffusion and migration of ions, and time-dependent behavior.

Description of the system
Here we describe variables that characterize the system's status.Consider a system with two phases: a solid phase s (electrode) and a liquid phase l (electrolyte), as shown in Figs.1a and b for potentiostatic and galvanostatic loading, respectively.The model is also valid for stress-free, single-phase solid electrolytes without grain boundaries and single ion conduction.Consider the following general reaction at the electrode-electrolyte interface where O i and R i are general oxidized and reduced states, respectively, and s o,i and s r,i are corresponding stoichiometric coefficients.The electrode-electrolyte interface will evolve due to the reaction, forming a potentially complex morphology.The location of the interface is described by an implicit function called the phase field variable ξ.The phase field variable changes from 1 in the bulk electrode to 0 in the bulk electrolyte and has a smooth transition across the interface, as shown in Fig. 1c.There is an intrinsic length scale ℓ associated with the diffuse interface.As shown in Fig. 1d, the electric potential ϕ typically shows a jump at the interface due to the electrical double layer.Here, we focus on systems with a size much larger than the Debye length; therefore, we do not explicitly resolve the double-layer structure.Instead, we define two electrostatic potential fields, ϕ s , and ϕ l , to implicitly capture the potential jump.As shown in Fig. 1d, both ϕ s and ϕ l are defined everywhere in the system but have the same value as ϕ in the corresponding bulk phases.Following the definition in the Kim-Kim-Suzuki or the grand potential models [21,23], they are called virtual electric potentials, which are only physical in the corresponding phases.
Considering an N -component system, the concentration of species i in phase α is denoted as c α i .We assume the two phases have the same molar volume V m ; in other words, N i=1 c α i = 1/V m , ∀α.This means only N − 1 concentrations are independent.By choosing a chargeneutral specie as the N -th component (charge number z N = 0), we define the diffusional electrochemical potential as μα i = μα i − μα N , where μα i = µ α i + F z i ϕ α is the electrochemical potential, µ α i is the chemical potential, z i is the charge number and F is the Faraday constant.Concentrations c s i and c l i are called virtual concentrations [21] (see Fig. 1e).The total concentration c i , is a combination of c s i and c l i by a mixture rule c i = p(ξ)c s i + (1 − p(ξ))c l i , where p(ξ) is an interpolation function given below.Generally, there is only one independent concentration [21,24], so we extend the classical local equilibrium constraint, which requires equal diffusional potential, to the electrochemical system as equal diffusional electrochemical potential μi = μs i = μl i to eliminate one extra degree of freedom.However, this local equilibrium constraint must be fulfilled at every grid point in the system, leading to higher computational costs.To solve this problem, we adopt a grand potential formulation [23,8] with the diffusional electrochemical potentials μi as independent state variables.In this way, the local equilibrium condition is fulfilled naturally.The local equilibrium constraint is used to construct the phase field model and does not mean the interface is at local equilibrium.It should be noted that we use the grand potential formulation to introduce the local equilibrium constraint and simplify the theoretical derivation.However, it is the Helmholtz free energy of the system that is minimized.A total electrostatic potential could be deduced from the equal diffusional electrochemical potential condition and the mixture rule for the total concentration; however, an explicit mixture rule for the total electrostatic potential is unnecessary.In summary, the status of the system can be uniquely described by the phase field ξ, the concentrations {c i } N −1 i=1 (natural variable for the Helmholtz free energy) or the diffusional electrochemical potential { μi } N −1 i=1 (natural variable for the grand potential), and the electrostatic potentials ϕ α , α = s, l.

Basic phase field theory
The electrical Helmholtz free energy of the system is where the grand potential of the system is defined as where where ω α , ⃗ D α and ⃗ E α are the grand potential density (J m −3 ), the electrical displacement (C m −2 ), and the electric field (V m −1 ) in phase α, respectively.Note that in previous models [5,12,10,8,14,17], the last term in Eq. 4 is neglected so ω α e = ω α .Here we keep the general form for consistency.Note here we use a grand potential formulation to decouple the bulk free energy from the surface energy; therefore, it enables the usage of a larger diffuse interface width [23,8].
The second law of thermodynamics requires (see Appendix A for derivation) where the total concentration (mol m −3 ) c i = −δΩ/δ μi fulfills the mass conservation equation where ⃗ J i is the total mass flux (mol m −2 s −1 ) of component i, which is defined by certain mixture rule from the phase mass flux ⃗ J α i .The charge densities are ρ α = ∂ω α /∂ϕ α .To fulfill the second law of thermodynamics (Eq.5), we can have the following constitutive relations from non-equilibrium thermodynamics [25] Here we define the variational overpotential [9,8] where n is the number of charges transferred (see Eq. 1).The constant before j(η) in Eq. 7 merely makes j the interfacial current density (A m −2 ).Here j(η) describes the general reaction in Eq. 1 and can be determined theoretically or experimentally, as long as it fulfills the second law (Eq.5 and using Eqs.7 and 10): Actually, Eq. 11 is fulfilled by typical choices of j(η): linear kinetics j(η) = kη [2,3,4,5], where k is a positive constant, or the Butler-Volmer-like kinetics [8] where j 0 is the exchange current density, and α a and α c are the anodic and cathodic charge transfer coefficients, respectively.In Eq. 9, we assume the electric field responds to any changes in the system instantaneously.This is valid since the characteristic timescale of the response of the electric field is proportional to the reciprocal of the speed of light, which is much smaller than any characteristic timescale of the interface movement or diffusion.We recover Gauss's law in Eq. 9.
For a system size much larger than the Debye length, this equation is typically replaced by the charge neutrality equation ρ α = 0 [26] and the charge conservation equation.Actually, the charge conservation equation can be derived by taking the time derivative of Eq. 9 and using Ampére's circuital law [27] (neglecting magnetic effects): where ⃗ i α is the current density (A m −3 ).Note that this treatment still obeys the second law in Eq. 5.
The kinetic equation Eq. 7, the diffusion equations Eq. 6 and Eq. 8 (contain both diffusion and migration of ions), and the charge conservation equations Eq. 13 construct the fundamental equations in our theory.It is worth noting that until this point, no material-specific assumptions, except the energy form given in Eq. 3, are made.This makes our theory very general and physically rigorous.In the next section, we will consider a particular case of binary electrolyte.

Binary electrolyte
Assume a binary salt dissolves in the solvent M ν + X ν − ↔ ν + M z + + ν − X z − .Four species are in the system: cation +, anion −, solvent N , and electron e − .Note we do not need to define the metal as a separate species since it can be viewed as a combination of a cation and an electron.The solvent component N is chosen to define the diffusional electrochemical potential.We assume electronic conduction dominates in the electrode while ionic conduction dominates in the electrolyte.In other words, we treat the electrode as a primarily electric conductor and the electrolyte as a primarily ionic conductor: We assume a simple interfacial reaction: M z + (l) + ne − (s) ↔ M (s), where n is the number of charges transferred per reaction (n = z + due to charge conservation).The current density in the electrode equals the electric current density ⃗ i s = ⃗ i s e = −σ s e ⃗ ∇ϕ s , and the current density in the electrolyte equals the ionic current density For simplicity, we made the following assumptions: (1) the off-diagonal elements of the mobility matrix are zero M ij = M i δ ij (no summation over i, δ ij is the Kronecker delta); (2) the ionic diffusivity in the electrode is much smaller than in the electrolyte D s + ≪ D l + ; (3) the concentration in the electrode c s + is a constant.With these assumptions, the mass conservation equation (Eq.6 and Eq.8) and the ionic charge conservation equation (Eq.13, α = l) can be derived as (see Appendix B for details) where − is the (electrolyte) ambipolar/effective diffusivity, and the ionic current in the electrolyte is The electric charge conservation equation (Eq.13, α = s) is In summary, Eqs. 7, 14, 15 and 17 are the governing equations.Until now, we haven't assumed any form of the Helmholtz free energy density or the relation between μi and c α i .In the next section, We will derive an expression for the overpotential using the free energy for the ideal solution.

Ideal solution free energy
To proceed, we need a constitutive relation of the free energy density.For simplicity, we assume an ideal solution of the electrical Helmholtz free energy densities (J m −3 ) where x α + = V m c α + is the molar fraction, A α and B α are materials parameters that can be determined experimentally.The grand potential densities ω α can be derived from Eq. 18 [23].With this, the variational overpotential in Eq. 10 can be derived as where 4 since this term is typically orders of magnitude smaller than ω α .Note that we recover the classical definition of the overpotential [28].The derivative of the interpolation function, p ′ (ξ), can be viewed as a regularized delta function, which merely restricts the quantity to the interface.The surface energy term in Eq. 19 is the diffuse interface representation of the capillary overpotential.

Relation to classical electrochemical quantities
In the classical electrochemical theory [29,28], the overpotential is only defined at the interface.It can be related to the variational overpotential in Eqs. 10 and 19 by integrating along a coordinate s perpendicular to the interface Similarly, the classical interfacial current density (A m −2 ) can be related to the corresponding phase field interfacial current density (A m −2 ): The total applied current (A) can be calculated as where the integral domain is the whole system.The capacity (Ah) can be calculated as the time integral of I: I dt.

Practical aspects
Higher-order kinetic correction For highly nonlinear kinetics or large overpotential, the phase field evolution equation in Eq. 7 can deviate from the desired sharp interface result j classical = j(η classical ).To improve the accuracy of the phase field kinetics, we expand j(η) into higher-order terms and match it with the sharp interface result to introduce the higher-order kinetic correction (see Appendix C for derivation) where η a equals the terms in the second bracket in Eq. 19.Note that this kinetic correction is not needed for linear kinetics since all higher-order terms are zero.The higher-order terms are also zero for the particular choice of the interpolation function p(ξ) = ξ.However, linear interpolation is generally not preferred as the equilibrium can shift from ξ = 0 and ξ = 1 [30], the overpotential is not localized near the interface (see Eq. 19), and it is not compatible with the driving force extension given below.In practice, we must truncate the expansion at a finite number N k .The choice of N k depends on the form of the reaction kinetics and the maximum overpotential.Details will be discussed below.Note that all the even order corrections are zero for Butler-Volmer kinetics with α a = α c = 1/2 and Marcus kinetics.In addition, a stable phase field profile can benefit from this correction at a large overpotential.

Galvanostatic condition
The galvanostatic condition (Fig. 1b) can be achieved by adjusting the applied voltage ϕ s = V app with a bisection algorithm [31] to match the current calculated from Eq. 22 to the applied current I app .If the spatial variation of the electric potential of the electrode ϕ s (Eq.17) is needed and the conductivity is constant, we can impose a constant charge flux through the Neumann boundary condition to achieve the galvanostatic condition.

General load profiles
The time-dependent loading, either the applied voltage V app (t) or the applied current I app (t), can be considered by a temporal discretization t k , k = 0, 1, • • • .For each timestep t k or stage (if a multi-stage time stepper is used), the problem can be treated the same as the case with a fixed applied voltage V app (t k ) or applied current I app (t k ).

Driving force extension
The grand potential formulation enables a thick interface width; however, a thick interface is not guaranteed to be numerically stable.The driving force, η a , is typically orders of magnitude larger than the surface energy term in Eq. 19.This leads to an unstable phase field profile and restricts the grid size to nanometers or even sub-nanometers.To enable a micrometer grid size, a driving force extension method is used to modify the overpotential in Eq. 19 as where P is the driving force extension operation that maps the driving force η a to a constant perpendicular to the interface.Details on the driving force extension method can be found in [22].The projection P keeps the first-order asymptotics of η a .Strictly speaking, this numerical treatment is non-variational.

The limiting current
The ideal solution free energy given in Eq. 18 requires a positive concentration.Numerically solving the equations with an applied charging current density close to the limiting current (typically > 97% of j lim ), the virtual concentrations can be negative in the nonphysical region (c l i in the electrode).To solve the numerical issue caused by a negative virtual concentration, we introduce a regularization whenever ln x α i is calculated: ln x α i → ln(max(ε, x α i )), where ε ≪ 1 is a small positive number.Note that even if the exponential function in the Butler-Volmer equation is combined with the logarithmic function in the overpotential, there are still terms with a negative power of x α i , forbidding a negative virtual concentration.Another way to handle the limiting current is to correct the Butler-Volmer equation as given in [13].In this work, we use the simple regularization given above.
It should be noted that charge separation needs to be considered [32,33] to correctly account for the limiting current and over-limiting current behavior.In this case, the charge neutrality assumption employed in this work cannot be used.Gauss's law in Eq. 9 should be used instead of the charge neutrality equation to resolve the electrical double layer explicitly.This will be addressed in a future publication.

Results
In this section, we first show the effect of the driving force extension method given above.Then, we provide a comprehensive verification of the proposed general theory, including (a) the equilibrium behavior, (b) the capillary effect, (c) the kinetic reaction, (d) the coupling between diffusion and migration of ions, and (e) the time-dependent behavior.We apply the proposed model to the electrodeposition of lithium metal.The electrolyte is 1M lithium hexafluorophosphate (LiPF 6 ) in ethylene carbonate (EC).We assume the molar fraction of Li + in the electrode is x s + = 0.99.The temperature is T = 298.15K. Materials parameters are: charge numbers z + = 1, z − = −1, the number of charges transferred n = 1, the molar volume V m = 1.302 × 10 −5 m 3 mol −1 [34], the equilibrium potential of Li|Li + is ∆ϕ eq = −1.26V [35], from which the free energy related parameters are determined as A s = 4.5920RT , A l = 0, B s = 44.4488RT, B l = 4.3282RT (R is the gas constant), the Li + diffusivity in electrolyte D l + = 3.2 × 10 −10 m 2 s −1 [36], the cation transference number t + = 0.5, and the surface energy σ = 1.176J m −2 [37].For Butler-Volmer kinetics, the exchange current density is j 0 = 30 A m −2 [37], and the symmetric coefficient is α = 0.5  The region above is unstable.With the driving force extension, there is no unstable region.The circular point shows an example.[38].For Marcus and Marcus-Hush-Chidsey kinetics, the exchange current density is j 0 = 104 A m −2 , and the reorganization energy is λ = 10RT [39].Unless explicitly mentioned, Buter-Volmer kinetics is used.The bulk electrolyte concentration is c 0 = 1 m.The diffuse interface width is ℓ = 1.5∆x.In this work, we focus on a half-cell and choose the electric potential on the right side of the electrolyte (see Fig. 1a) as the reference ϕ l = 0 V.The electrode is assumed to be equipotential, whose electric potential equals the applied voltage ϕ s = V app .The half-cell voltage is then E cell = ϕ l − ϕ s = −V app .Numerically, we use the finite difference method for spatial discretization with second-order central differences for the derivatives.Forward-Euler temporal discretization is used for the phase field equation Eq. 7, and backward-Euler is used for the concentration equation Eq. 14.The resulting linear systems and the Poisson equations are solved by the Bi-CGSTAB iterative solver [40] with a multigrid preconditioner.

Driving force extension
The surface energy term in Eq. 19 can stabilize the phase field profile during evolution.Without the driving force extension in Eq. 24, the term related to η a typically has a destabilizing effect.To keep a stable traveling wave profile, it is required that these two terms have a similar order of magnitude.The magnitude of the surface energy term is approximately 6V m σ/(nF ℓ).Considering a practical overpotential, η a = 200 mV, the maximal interface width is restricted to be ℓ ∼ 6V m σ/(nF η a ) = 4.8 nm (this is just a numerical requirement, not a physical one).The grid size will be even smaller to properly resolve the diffuse interface.The maximum applied voltage with a stable phase field profile for a given grid size is shown in Fig. 2. For practical applied voltages, the stable grid size is on the order of nanometers.This is a significant limitation for simulating practical system sizes that range from micrometers to millimeters.We apply the driving force extension method given above to remove this limitation.As shown in Fig. 2, for a practical load of 200 mV, the phase field profile is only stable with a 1 nm grid size without the driving force extension while stable with a 1 µm grid size with the driving force extension.Note that three orders of magnitude increase in grid spacing results in six orders of magnitude increase in timestep size for explicit time stepping.This leads to a significant improvement in computational speed, estimated to be 12 and 15 orders of magnitude for 2D and 3D, respectively, for explicit

Equilibrium behavior
We first verify if the proposed model can capture the correct equilibrium behavior.The equilibrium concentration and the cell voltage are related by the Nernst equation [28] where E • − cell = −∆ϕ eq is the equilibrium cell voltage at c l + = c • − = 1 m.We use step voltage testing on a system given in Fig. 1a.The system size is 100 µm, and the grid size is 1 µm.We solve the proposed model in 1D (planar interface) with zero-Neumann boundary conditions for both the phase field and the concentration.For the electrostatic potential, the Dirichlet boundary condition is used on the right ϕ l = 0 V and the zero-Neumann boundary condition on the left.Initially, the interface is located in the middle of the system, the electrolyte has a homogeneous concentration c l + = 1 m, and the applied voltage is We then step the applied voltage with a size of −5 mV, as shown in Fig. 3a.The interface moves, and the concentration relaxes toward a new equilibrium value determined by the applied voltage; see Fig. 3b.After complete relaxation (here defined to be the changes in average concentration c l + being less than 7.68 × 10 −8 m), we record the average c l + and the applied voltage V app and repeat the voltage stepping.The recorded concentrations and applied voltages are plotted in Fig. 3c and compared with the Nernst equation (solid line).This verifies that our model correctly captures the equilibrium behavior.

Capillary effect
Considering an infinitely long cylindrical electrode with radius r, the equilibrium concentration and the cell voltage are related by a modified Nernst equation, accounting for the capillary effect To model this problem, we solve the proposed model in 1D cylindrical coordinates.To exaggerate the capillary effect, we model a small cylindrical system of 200 nm radius with the electrode located at the center and the electrolyte elsewhere.We apply a voltage The equilibrium concentration for a planar interface (r → ∞) will be 1 m.Due to the capillary effect, the equilibrium concentration with the cylindrical electrode will deviate from 1 m.Boundary conditions are the same as in the Equilibrium behavior section.The grid size is ∆x = 0.1 nm.The electrode radius varies from 10 nm to 150 nm, and the average concentration is measured after full relaxation.A stationary interface is used to fix the electrode radius.This is achieved by freezing the phase field ξ.As shown in Fig. 4, the phase field results compare well with the analytical solution.Moreover, as the electrode radius increases, the equilibrium concentration approaches the equilibrium concentration for a planar interface (shown by the gray dashed line).Note that after full relaxation, the interface will remain stationary even if we let ξ evolve.This verifies that our model captures the correct capillary behavior, critical in many applications such as Li dendrite growth.
To verify the results for a broad overpotential range, we manually fix c l + = 1 m and do not evolve the phase field ξ (can be viewed as in a moving reference frame).We model a planar interface and apply linear potential sweep with a rate of −1 mV s −1 .The settings are the same as in the Equilibrium behavior section.Higher-order kinetic correction given in Eq. 23 is used but truncated at N k .During the voltage sweeping, the overpotential and current density are calculated from Eqs. 20 and 21, respectively, and are shown in Fig. 5 for both results without the kinetic correction (N k = 0) and with a sufficient order of kinetic correction.
For the Butler-Volmer equation, the result without correction (Eq.7) and with higherorder correction (Eq.23, truncated at N k = 11) are compared in Fig. 5a.It can be seen that the kinetic correction improves the accuracy for large overpotentials.Figure 5b shows the relative error of the current density at two overpotential values as a function of the kinetic correction order N k .The error reduces quickly as more orders of corrections are used.For smaller overpotential, fewer corrections are needed.For the rest of the test cases in this paper, Butler-Volmer kinetics with N k = 11 is used.The results for the Marcus and the Marcus-Hush-Chidsey models are shown in Figs.5c and d, respectively.The behavior is similar to the Butler-Volmer case.However, the Marcus kinetics is more nonlinear, so more correction terms are needed.For the Marcus-Hush-Chidsey model, we approximate the Marcus-Hush-Chidsey model by a 13-order polynomial (see Appendix D), so N k = 13 gives the exact result.This verifies that our model can capture general reaction kinetics correctly, even at high overpotentials.

Coupling between diffusion and migration
To test the coupling between diffusion and migration of ions, we verify our model against the analytical result using a sharp interface model where boundary conditions are applied at the interface, at quasi-steady-state (∂c l + /∂t = 0) for a planar interface.The geometry can be seen in Fig. 6   coupled way.To be consistent with the analytical solution, the interface is assumed to be stationary (not evolving ξ).The phase field results for the electrolyte concentration c l + and the electrostatic potential ϕ l are compared with the analytical solution [37] in Figs.6a and  6b, respectively.Note that if a single electrostatic potential field ϕ is used, there will be an electrostatic potential jump across the electrode-electrolyte interface from ϕ l ∼ −30 mV to ϕ s = −1460 mV, which may be difficult to be resolved accurately in a numerical solver.For higher dimensions, the potential jump will differ from point to point on the interface, which traditional methods overlook.With the two-electrostatic-potential model proposed here, the jump can be correctly handled without explicitly resolving the potential jump.This verifies that the proposed model gives the correct coupling behavior between diffusion and migration of ions in the electrolyte.

Transient behavior
To examine the time-dependent behavior of the proposed phase field model, we consider the initial transient of a system with a flat electrode-electrolyte interface under galvanostatic conditions.Initially, the system is under open circuit condition, and the electrolyte is homogeneous with a concentration c 0 = 1 m.The geometry is shown in Fig. 7. Boundary conditions are the same as in the previous section.At t = 0, we apply a constant charging current of 95% of the limiting current j lim = 2F D l + c 0 /L = 617.5A m −2 , where L = 100 µm is the system size.With time, the cation depletes close to the interface and eventually reaches a quasi-steady state.The simulated Li + concentration profile is compared with the analytical solution for the sharp interface model [37] at various times in Fig. 7.This indicates that our model can quantitatively predict time-dependent behaviors.

Discussion
In summary, the proposed phase field method can successfully model various aspects of electrochemical systems with realistic materials parameters and practical system sizes.The model is fully variational and computationally tractable.Here we discuss possible extensions of the proposed phase field model.

Surface energy anisotropy
Surface energy anisotropy can be included in the gradient coefficient κ in Eq. 3 by letting κ depend on the orientation of the interface through the normal to a level-curve of ξ, n = − ⃗ ∇ξ/| ⃗ ∇ξ|.The variational overpotential (Eq.10) then needs to include extra terms related to the derivative of κ to ξ.This is straightforward and does not change the applicability of the proposed theory and is, therefore, not elaborated here.Details on introducing anisotropy can be found in e.g., [42,43].For strong anisotropy, convexification can be employed [44] or the curvature dependence of the interfacial energy can be introduced [45,46,47,48].

General free energy forms of the electrolyte solution
The proposed theory is formulated to handle general forms of free energy density.The ideal solution given in Eq. 18 is used here for simplicity.Consider a general form of the electrochemical potential for concentrated solutions, such as where a i is the activity that depends on all the independent concentrations c i .The corresponding free energy density can be determined as f ({c i } N −1 i=1 ) = N i=1 c i μi , from which the grand potential density ω can be determined.The free energy density can also be defined by the Redlich-Kister expansion [17].The derivation of the governing equations and variational overpotential can then proceed similarly.

Field-dependent materials parameters
Materials parameters can depend on certain field variables.For example, the ionic diffusivity may depend on the local temperature and concentration.This may introduce additional numerical difficulties, but the proposed theory can be applied to these cases without any change.

Other physical phenomena
One of the main advantages of the phase field method is that it is easy to add new physics to the model.For example, to account for the stresses in a solid electrolyte, where the strain energy can be included in the total free energy in Eq. 3, the relevant evolution equations can be derived following the same procedure.If the evolution of the thermal field is needed, the entropy formulation [30,49,50] that maximizes the total entropy can be used.To investigate polycrystalline metal anode or the grain structure of solid electrolytes, phase field models with multiple phases [24] can be used.

SEI effect
SEI can form at the electrode-electrolyte interface, especially for alkali metal anodes.Since the SEI thickness is typically on the nanometer scale, explicitly considering the SEI as a separate phase limits the simulated system size e.g. in [51].However, the effect of the SEI can be incorporated into the proposed model through parameter modification.For example, the lower ionic diffusivity of the SEI can be included by making the ionic diffusivity dependent on the order parameter [15].The stress due to a SEI layer can be incorporated as surface stress.Implicit coupling of the SEI effect into our model will be addressed in a future publication.

Conclusion
In summary, we have introduced a general theory that can be used for quantitative phase field modeling of electrochemical systems that employ realistic materials parameters and practical system sizes.The model is fully variational, guaranteeing that the system's energy must decay during evolution in an isothermal system.The classical definition of the overpotential is recovered for a binary electrolyte and ideal solution free energy.We employ a grand potential formulation to decouple the bulk free energy from the interfacial energy, enabling a much larger diffuse interface width.The classical local equilibrium constraint is extended to the electrochemical case using the diffusional electrochemical potential.In addition, a driving force extension method is applied to increase the grid size by three orders of magnitude and the explicit timestep size by six orders of magnitude compared to standard methods that can be used in a phase field calculation.Thus, the proposed phase field model can be used to efficiently simulate electrochemical systems of practical size with realistic materials parameters.
By incorporating the reaction kinetics in the evolution equation for the phase field parameter, the proposed model can handle general reaction kinetics without additional model parameter calibration.This treatment allows a correct treatment of the capillary overpotential.A higher-order kinetic correction is proposed to accurately handle a range of reaction kinetics models for any overpotential, such as Butler-Volmer kinetics, Marcus kinetics, and Marcus-Hush-Chidsey kinetics.
Both the electrostatic potentials in the electrode and the electrolyte are chosen as independent variables.Compared with previous models where only one electrostatic potential is used, this model allows for more freedom in capturing the correct potential jump at the electrode-electrolyte interface, including a variation of the jump along the interface, without explicitly resolving the potential jump.
The model has been verified using several tests: • The model gives the proper equilibrium behavior described by the Nernst equation.
• The critical role of interfacial energy in the overpotential is captured.
• General reaction kinetics at the electrode-electrolyte interface are captured.
• The coupling between diffusion and migration of ions as given by classical theory is recovered.
• The time evolution of the concentrations in the electrolyte agrees with analytical models.
The proposed phase field framework is very general.It allows applications in a broad range of materials systems, e.g., general electrolyte solution (multicomponent, general charge number), general reaction kinetics, general free energy form for the electrolyte solution (concentrated solutions), and general load profiles (potentiostat, galvanostatic, or timedependent loading).The nature of the phase field model allows easy extension to include additional physics, like mechanical and thermal effects.

Figure 1 .
Figure 1.Schematic illustration of the electrochemical system for potentiostatic (a) and galvanostatic (b)conditions.E cell = −V app is the half-cell voltage.Illustration of field variables to describe the status of the system: (c) the phase field variable ξ varying across a length scale called the diffuse interface width l, (d) the electrostatic potential ϕ and the virtual electrostatic potentials ϕ s and ϕ l for each phase, and (e) the concentration c i of species i and the virtual concentrations c s i and c l i for each phase.
an interpolation function, κ = 6σℓ and m = 3σ/ℓ are model parameters related to the surface energy σ and the diffuse interface width ℓ, and[25,11] where σ s e is the electric conductivity of the electrode.The electrical charge densities in each phase are ρ s e = −F c s e − = −F z + c s + and ρ l e = 0.The ionic charge densities are ρ s ion

Figure 2 .
Figure 2. Driving force extension.The solid line shows the maximum applied voltage (versus open-circuit voltage) that gives a stable phase field profile for a given grid size without the driving force extension (DFE).The region above is unstable.With the driving force extension, there is no unstable region.The circular point shows an example.

Figure 3 .
Figure 3. Equilibrium behavior.(a) The applied voltage step history.(b) The Li + concentration relaxation history compared to the equilibrium concentration determined by the Nernst equation.(c) The equilibrium concentration as a function of the applied voltage (versus open-circuit voltage) compared with the Nernst equation.

Figure 4 .
Figure 4. Capillary effect.The equilibrium electrolyte concentration c l + as a function of the electrode radius r for an applied potential equals the open-circuit voltage for a planar interface.The gray dashed line shows the equilibrium concentration for a planar interface (r → ∞).

Figure 5 .
Figure 5. Reaction kinetics.(a) Comparison of the calculated interfacial current density for the Butler-Volmer model without kinetic correction (N k = 0) and with kinetic correction truncated at N k = 11 with the analytical Butler-Volmer equation.(b) The effect of the kinetic correction order on the relative error in the current density at two overpotential values for the Butler-Volmer case.Comparison of the calculated current density with analytical expression with and without kinetic correction for Marcus (c) and Marcus-Hush-Chidsey (d) models.
. Dirichlet boundary conditions for the concentration c l + = c 0 and the electrostatic potential ϕ l = 0 are used for the right boundary.We apply a voltage V app = −E • − cell − 200 mV on the electrode.The rest of the settings are the same as the Equilibrium behavior section.Governing equations Eqs.14, 15 and 23 are solved in a fully

Figure 6 .
Figure 6.Concentration and potential response of a planar interface to an applied voltage at quasi-steadystate.Comparison of the phase field results with analytical solutions for (a) the concentration field c l + and (b) the electrostatic potential field ϕ l .Note that the virtual concentration c l+ has no physical meaning in the electrode, and it is merely a solution to chemical equilibrium.

Figure 7 .
Figure 7. Initial transient of an initially homogeneous electrolyte in response to an applied current of 95% of the limiting current.The phase field calculated electrolyte concentrations c l + (solid lines) are compared with the analytical solution (dashed lines) at various times.Note that the virtual concentration c l + has no physical meaning in the electrode.