Moreover, even if a numerical technique is employed, it requires an explicit tracking of the interface boundaries [81, 163].
This is a cumbersome theoretical challenge, especially since the interface evolves continuously and with complex geometrical shapes [164].
Finally, the responses in a local electric field due to the dynamic change in the inclusion shape are neglected [82, 83].
In the present work, a phase-field model is derived to investigate the temporal evolution of the inclusions propagating in the conductor.
The elegance of the phase-field method lies in its ability to simulate moving boundary problems without having to track the interfaces explicitly. In the following section, a free energy functional is described to track a morphologically evolving inclusion in a conductor.

§ FREE ENERGY FUNCTIONAL

[Schematic diagram of the simulation setup describes the inclusion in the matrix subjected to external electric potential]A schematic diagram of the simulation setup describes the inclusion with the order parameter $c=1$, in the interconnect material $c=0$, which is subjected to the electric potential $+\phi_\infty$$\phi_\infty$electrical potential boundary condition at the anode, and $-\phi_\infty$ at the cathode.

A phase-field model is employed to examine the electromigration-driven dynamics of transgranular inclusion. The schematic of an inclusion subjected to an external electric field is presented in Figure  <ref>. A conserved order parameter $c$ is introduced to demarcate the matrix denoted by $c = 0$ and the inclusion where it assumes a value of $c = 1$. The traditional sharp interface between the matrix and the inclusion is replaced by a narrow region of finite thickness. In this region, the order parameter varies smoothly between 0 and 1. In this way, the variable $c$ fulfills a dual role in tracking the species concentration and the interface between the matrix and the inclusion. The smooth interface provides a remarkable advantage in phase-field modeling to avoid the tedious task of explicit tracking as opposed to the sharp interface counterpart.

The free energy of the system $F(c,\boldsymbol{\nabla} c)$ is expressed as a function of the order parameter and its gradient as

\begin{equation}
F(c,\boldsymbol{\nabla} c) = \int_{V_\Omega} \left[ f(c) + \frac{\kappa}{2} \vert \boldsymbol{\nabla} c \vert^2 \right] \textrm{d}{V_\Omega}, \label{eq:EMfreeEnergyFunctional}
\end{equation}

where $f(c)$ represents the bulk free energy, $\kappa$ is the gradient energy coefficient, and $ \vert \bullet \vert$ denotes the norm of the vector.

It is worth noticing that the choice of the free energy function, $f(c)$, can produce a significant effect on the physical behavior of the interface. Hence, it should be selected appropriately. The regular solution model derived in section <ref> is a well-known free energy function for a broad range of applications. In that, the position of the two wells is decisive to obtain equilibrium values based on optimal energy state.
For symmetric energy curve ($X_1=0$) and $T=T_{\textrm{ref}}$ in Eq. (<ref>), the minima of the free energy are the equilibrium values [For asymmetric curves, a Maxwell construction can be employed to determine the equilibrium values, which is further discussed in details in section <ref>. Furthermore, these values are pertinent to flat interfaces. Alternatively, the curvature correction should be added according to Gibbs-Thomson law [165].]. For instance, at $X_2=2.6$, two equilibrium values are $c^L=0.124$ and $c^H=0.876$. In addition, $X_2=2$ represents the critical limit below which there is no concavity and the free energy plot resembles to Figure <ref>(c). As the $X_2$ increases above a critical limit, the minima starting to drift apart progressively. For a very high $X_2$, the two equilibrium values coincide $c^L=0$ and $c^H=1$. These values are convenient to indicate the bulk phases [139, 166, 167, 142].
However, the numerical calculations are unstable at these points due to logarithmic terms in Eq. (<ref>). An alternative form of the free energy can be expressed as Landau polynomials [139, 123] of the form,

\begin{equation}
f^L(c) = \sum_{l=0}^{N_l} L_l c^l,
\end{equation}

where $l$ indicates the index, $N_l$ represents the last term, and $L_l$ is the coefficient of the $l$-th term of the Landau polynomial. The simplest form of the polynomial is the double-well function, $f^{dw}(c)=c^2(1-c)^2$.

The double-well function represents an approximation of the Van der Waals [168, 169] near the critical point, and has been employed extensively in the phase-field models. However, when the model is developed solely for interface tracking purposes, this has led to the frequently observed spontaneous volume shrinkage phenomenon. Whereby, the high volume phase allows significant infiltration into the low volume phases and can ultimately cause the complete disappearance of the lower volume phases [170]. In addition, the movement of the lower volume phase inside the higher one may further enhance this effect. Therefore, to minimize these losses over the duration of a simulation, it requires a reconsideration of the free energy function. An alternative energy function for interface tracking applications in the form of the double-obstacle can be considered. The diffuse interface of the order parameter follows cosine function in the double-obstacle, while hyperbolic-tangent is expected in the former case, see Appendix <ref>. This provides an extremely controlled diffuse interface width and less dispersion of inclusion volume. Consequently, lower change in inclusion volume leads to a stoppage of the spontaneous volume shrinkage phenomena. In addition, this function is adequate to track larger movements of the lower volume phase. Therefore, the obstacle-type function may prove useful for interface tracking applications of the phase-field model where the nature of the simulated phenomena introduces phase continuity concerns like the inclusion migration in the metal conductors.

In addition, the obstacle-type free energy gains a computational advantage, by allowing the solution needs to be determined merely in the interface region, over the double-well type formulation [134]. The obstacle-type free energy density $f^{\textrm{ob}}(c)$$f^{\textrm{ob}}$obstacle-type free energy density, as shown in Figure <ref>, is considered with equal minima at $c=0$ and 1:

\begin{equation}
f(c)= f^{\textrm{ob}}(c) = X_A c (1-c) + I(c),
\end{equation}

where $X_A$$X_A$barrier height of the obstacle-type free energy density sets the barrier height of the free energy density, and $I(c)$$I$indicator function for the obstacle-type free energy density is the indicator function expressed as

\begin{eqnarray}
I_{[0,1]} =     \left\{ \begin{array}{ccl}
0, & \mbox{  for}& 0 \le c \le 1 , \\
\infty,  & \mbox{  for} & c<0  \mbox{   \&   } c>1.
\end{array}\right.%
\end{eqnarray}

[Double-obstacle type free energy density $f^{\textrm{ob}}(c)$ as a function of $c$]Double-obstacle type free energy $f^{\textrm{ob}}(c)$ as a function of $c$. The function is symmetric about center $c=0.5$ and has two free energy minima at 0 and 1.

The morphological evolution of inclusion can be recognized by the competition between the electromigration force and the capillary force. The capillarity alone prefers circular shape inclusion. While the external electric field instigates species transport, which leads to shape alterations. In addition, the capillary force seeks uniform curvature and consequently diffuses the species to reduce any curvature gradient along the inclusion surface. Therefore, the shape of the inclusion is governed by the relative strength of these two forces, which can be incorporated into the diffusional fluxes.

§ DIFFUSIONAL FLUXES AND TIME EVOLUTION

The evolution of the diffusing species follows a continuity equation

\begin{equation}\label{eq:EMFluxEvolution}
\frac{\partial c}{\partial t} = - \boldsymbol{\nabla} \cdot \boldsymbol{J}_i,
\end{equation}

where $t$ denotes time, $\boldsymbol{\nabla} \cdot (\bullet)$ represents the divergence of the vector and $\boldsymbol{J}_i$$\boldsymbol{J}_i$flux of the diffusing species $i$ is the net flux of the diffusing species $i$. The flux $\boldsymbol{J}_i$ is expressed as a linear combination of the flux and driving force within the framework of irreversible thermodynamics as [69, 171]

\begin{equation}\label{eq:EMFluxVsChemicalPotential}
\boldsymbol{J}_i = -M_{ii}\boldsymbol{\nabla}(\mu + eZ_i\phi) -M_{ie}eZ_w\boldsymbol{\nabla}\phi
\end{equation}

where $M_{ii}$$M_{ii}$mobility of species $i$ due to the interaction with the species $i$ and $M_{ie}$$M_{ie}$mobility of species $i$ due to the interaction with the electrons are the mobilities related to the diffusivity of the species $i$ and interaction between species $i$ and electron respectively. $\boldsymbol{\nabla} (\bullet)$ represents the gradient of a scalar, $\mu$ is the chemical potential, $Z_i$$Z_i$valence of the diffusing metal species $i$ denotes the valence of the diffusing species, $Z_w$$Z_w$momentum exchange effect between the electrons and the diffusing species denotes the momentum exchange effect between the electrons and the species, $e$ denotes the electric charge and $\phi$$\phi$local electrical potential is the electrical potential. The term $(\mu + eZ_i\phi)$ combined is referred as the electrochemical potential [172]. To perceive the driving forces arising due to the imposed electric field, Eq. (<ref>) can be rearranged to give

\begin{equation}\label{eq:EMFluxModified}
\boldsymbol{J}_i = -M_{ii}\boldsymbol{\nabla}\mu - M_{ii}eZ_i\boldsymbol{\nabla}\phi - M_{ie}eZ_w\boldsymbol{\nabla}\phi
\end{equation}

The second term in the previous equation is the result of the direct electrostatic force, while the third term is the electron wind force, which reflects the cross-effect arising from the interaction between the diffusing species and the conducting electrons. In conductors, metal species are shielded by the negative electrons so that the direct electrostatic force is much less than the wind force [69, 173]. The dominance of the wind force is further corroborated by experimental observation which suggests the movement of the diffusing species in the direction of the electron flow [174, 68]. Therefore, the second and third terms are combined through effective valence, $Z_s$. Hence considering the effective valence in Eq. (<ref>) and substituting in Eq. (<ref>), the modified Cahn-Hilliard equation can be expressed as,

\begin{equation}\label{eq:EMCahnHilliardModified}
\frac{\partial c}{\partial t} = \boldsymbol{\nabla} \cdot \left(M_{ii}\boldsymbol{\nabla}\left(\mu +eZ_s \phi \right)\right)
\end{equation}

The chemical potential $\mu$ can be obtained from the variational derivative of the system free energy,

\begin{equation}
\mu = \frac{\delta F}{\delta c}.
\end{equation}

The surface diffusion is dominant compared to the other forms of mass transport mechanisms at the operating conditions. Therefore, the surface effective charge $Z_s$ and the atomic mobility $M_{ii}$ are restricted at the surface by selecting a bi-quadratic form in $c$, expressed as,

\begin{equation}
M_{ii}(c)= D_s f^{\theta}(\theta) c^2(1-c)^2,\label{eq:EMmobilityFunction}
\end{equation}

where $D_s$ denotes the surface diffusion coefficient$D_s$surface diffusion coefficient. Here $f^{\theta}(\theta)$$f^{\theta}$anisotropy function is the anisotropy function, given by [116, 115],

\begin{equation}
f^{\theta}\left(\theta\right) =\frac{\left(1+A \, \textrm{cos}^2\left(m\left(\theta+ \varpi\right)\right)\right)}{\left(1 +A\right)}. \label{eq:anisotropy}
\end{equation}

In the above equation, $A$$A$strength of anisotropy represents the strength of anisotropy, which describes the superiority of the maximum value of diffusion compared to the minimum along the surface. In addition, $m$$m$grain symmetry parameter is the parameter related to the grain symmetry. Specifically, $2m$ denotes the number of crystallographic directions of fast diffusion sites in the plane of inclusion migration. Therefore, based on a total number of locations for fast diffusivity, $m=1,\,2,$ and 3 are characterized by twofold, fourfold, and sixfold symmetry respectively, as shown in Figure <ref>, while $m=0$ indicates isotropic inclusions. Furthermore, $\theta = \tan^{-1} \big( \frac{\partial c}{\partial y}/\frac{\partial c}{\partial x}\big)$$\theta$angle formed by the local tangent at the inclusion surface is the angle formed by the local tangent at the inclusion surface and $\varpi$$\varpi$misorientation angle is the misorientation angle formed in the clockwise direction by the position of the maximum diffusion site with the perpendicular to the external electric field.

[Schematic of anisotropy in surface diffusivity of inclusions]Schematic of anisotropy in surface diffusivity with (a) twofold, m=1, (b) fourfold, m=2, and (c) sixfold, m=3 symmetrical inclusions. The double-headed arrows define positions of maximum surface diffusivity. A misorientation angle, $\varpi$ for each case is defined with respect to the direction of the perpendicular to the external electric field.

§ ELECTRICAL POTENTIAL

The flux of the charge carrier, i.e. electron, $\mathbf{J}_{e}$$\boldsymbol{J}_{e}$flux of the charge carrier is given by

\begin{equation}
\mathbf{J}_e = -M_{ei} \boldsymbol{\nabla}(\mu + eZ_i\phi) - M_{ee}eZ_w\boldsymbol{\nabla}\phi
\end{equation}

where $M_{ei}$$M_{ei}$mobility of electrons due to interaction with species $i$ and $M_{ee}$$M_{ee}$mobility of electrons due to interaction with the electrons are the phenomenological coefficients and $e$$e$electron charge is the electron charge. The current inside the conductor is entirely due to the imposed electron wind and the cross effect due to mass flux i.e. the first term is negligible. In addition, the timescale of relaxation of charge is much faster compared to the diffusion process. Hence, the current continuity equation translates into Laplace equation as

\begin{equation}
\boldsymbol{\nabla} \cdot \mathbf{J}_e = \boldsymbol{\nabla} \cdot [M_{ee}eZ_w\boldsymbol{\nabla}\phi] = 0.
\end{equation}

Comparing the previous equation with Ohm's Law, the conductivity function can be written as, $\sigma(c) = M_{ee}eZ_w.$ Therefore, the potential field is calculated from the Laplace equation as

\begin{equation}
\boldsymbol{\nabla} \cdot [\sigma(c)\boldsymbol{\nabla}\phi] = 0, \label{eq:EMLaplaceeq}
\end{equation}

where $\sigma(c)$$\sigma$electrical conductivity dependent on order parameter $c$ is the electrical conductivity dependent on order parameter $c$, interpolated between the inclusion $\sigma_{\textrm{icl}}$$\sigma_{\textrm{icl}}$conductivity of an inclusion and the matrix $\sigma_{\textrm{mat}}$$\sigma_{\textrm{mat}}$conductivity of the matrix as,

\begin{equation}
\sigma(c) = \sigma_{\textrm{icl}}h(c) + \sigma_{\textrm{mat}}[1-h(c)].
\end{equation}

Any smooth function that satisfies $h(c)\vert_{c=0}$ = 0 and $h(c)\vert_{c=1}$ =1 can be employed. Some functions that satisfy these properties are $h(c) = c, h(c) = c^2(3-2c),$ and $h(c) = c^3(10-15c+6c^2)$. The linear interpolation function ($h(c) = c$) is considered in the present work for computational convenience and the form of interpolation functions does not alter the results as long as the interface width is small.

§ RESEMBLANCE TO SHARP-INTERFACE DESCRIPTION

To facilitate comparison and physical significance, the relations of phase-field model parameters should correspond to sharp-interface theories. In this section, a linkage between the two methods is established.

§.§ Interfacial properties

In the phase-field model, the physical properties of the phenomenon such as interfacial energy $\gamma_s$$\gamma_s$interfacial energy and interfacial width  $\delta_s$$\delta_s$interface width are related to the double-obstacle barrier height $X_A$ and the gradient-energy coefficient $\kappa$.

The spatial gradient of the order parameter along the interface can be equated to the free energy density function to derive interfacial width, rearranging Euler-Lagrange Eq. (<ref>), of the form

\begin{equation}
{\textrm{d}x} = \sqrt{\frac{\kappa}{2f(c)}}{\textrm{d}c}, \label{eq:EMdcbydxAlongInterface}
\end{equation}

Substituting $f(c)$ from Eq. (<ref>) into Eq. (<ref>) and integrating along the interface from $c=0$ to $c=1$ $\Rightarrow$ $x=0$ to $x= \delta_s$. The interface width can be obtained as,

\begin{equation}
\delta_s = \pi \sqrt{\frac{\kappa}{2X_A}}. \label{eq:EMinterfaceWidth}
\end{equation}

$\pi$Archimedes' constantIt is important to note that this relation is applicable to obstacle type free energy density, where the interface width can be strictly related to $\kappa$ and $X_A$. However, for well-type free energy density, the diffuse interface reaches infinity. Therefore, to determine the interface width, an approximation should be considered by limiting the change in order parameter upto a few orders.

The system free energy expressed in Eq. (<ref>) identically vanishes in the bulk phases where $c=0$ and 1, while it has finite value at the interface. In other words, since the energy of the system is solely dictated by the movement of the interface, the interfacial energy can be expressed of the form,

\begin{equation}
\gamma_s = \int_{V_\Omega} \Big\{ f^{\textrm{ob}}(c) + \frac{\kappa}{2} \vert \boldsymbol{\nabla} c \vert^2 \Big\} \textrm{d}{V_\Omega}.
\end{equation}

This expression can be calculated numerically to obtain the interfacial energy for a complex system of the free energy densities [54]. Alternatively, the derivation of an analytical relation is also plausible for the simple systems considered for the presented case. To derive the expression, assuming that the spatial dependency of the order parameter $c$ is restricted to X-direction, writing one-dimensional form of the previous equation,

\begin{equation}
\gamma_s = \int_{-\infty}^{ \infty} \Bigg\{ f^{\textrm{ob}}(c) + \frac{\kappa}{2} \left( \frac{\textrm{d} c}{\textrm{d}x} \right)^2 \Bigg\} \textrm{d}x.
\end{equation}

Substituting the Euler-Lagrange Eq. (<ref>) into the above equation,

\begin{equation}
\gamma_s = \sqrt{2\kappa} \int_{0}^{1}  \sqrt{f^{\textrm{ob}}(c)}  \textrm{d}c. \label{eq:EMinterfacialEnergyINtermediate}
\end{equation}

Substituting $f^{\textrm{ob}}(c)$ from Eq. (<ref>) into Eq. (<ref>) and integrating, the final form of the interfacial energy can be expressed as,

\begin{equation}
\gamma_s = \frac{\pi}{4\sqrt{2}} \sqrt{\kappa X_A}. \label{eq:EMinterfacialEnergy}
\end{equation}

It is evident from these relations (<ref>) and (<ref>) that the physical properties of the system, the interface width and the interfacial energy can be recovered in a phase-field model. Furthermore, by manipulating the energy barrier height $X_A = 4 \gamma_s/\delta_s$ and the gradient coefficient $\kappa = 8 \gamma_s \delta_s/\pi^2$ desired values of interface width $\delta_s$ and the interfacial energy $\gamma_s$ can be achieved.

§.§ Determination of Mullins' constant

[Evolution of the sinusoidal amplitude with time.]Evolution of the sinusoidal amplitude with time. $A_0$ denotes the initial amplitude of a perturbation with frequency $k$. At time $t$, $A_t$ denotes the temporal-dependent amplitude of the perturbations. The points are obtained from phase-field simulations, which are approximated by a linear line. The slope of the line indicates $-D_s \delta_s \Omega \gamma_s k^4/(k_BT)$.

[Decay in the amplitude of the sinusoidal wave with time]Decay in the amplitude of the sinusoidal wave with time sequence from (a) to (c).

Surface diffusivity related to the Mullins' constant $B$Mullins' constant$\Omega$atomic volume$B (=D_s \delta_s \Omega \gamma_s/k_B T)$ is obtained by considering the dampening of a sinusoidal perturbation $ A_0$ sin($kx$) under surface diffusion [175], as shown in Figure <ref>. Here, $A_0$$A_0$amplitude of the perturbation at time $t=0$ represents the amplitude of the perturbation at time $t=0$, and $k$$k$frequency of the sinusoidal perturbation denotes the frequency of the perturbation.
According to the Mullins' law for the curvature-driven surface diffusion [176], the progressive dampening of the perturbation is expressed by

\begin{equation}
\frac{\partial y}{\partial t} = -\left( \frac{D_s \delta_s \Omega \gamma_s}{k_BT} \right) \frac{\partial^4 y}{\partial x^4}.
\end{equation}

The analytical solution of this equation, which is valid for $A_0 k<1$, provides the expression for the amplitude of the perturbation:$A_t$amplitude of the perturbation at time $t$

\begin{equation}
A_t(t) = A_0 exp \left[ \left( -\frac{D_s \delta_s \Omega \gamma_s}{k_BT} \right) k^4 (t) \right].
\end{equation}

Thus, the value of the parameter $D_s/k_BT$ is obtained by the plot of $\textrm{ln}(A_t/A_0)$ vs. $t-t_0$ as shown in Figure <ref>, which provides the slope = $\left( -\frac{D_s \delta_s \Omega \gamma_s}{k_BT} \right) k^4$.

§ NUMERICAL PROCEDURES

Eqs. (<ref>) and (<ref>) constitute the coupled partial differential equations to describe the inclusion dynamics under electromigration. It has been reported in Ref. [156] via formal asymptotic analysis, that at the sharp-interface limit i.e. when the interface width tends to zero, these coupled equations recover the motion by surface laplacian of the curvature and the electric potential. In addition, these equations are normalized for the computational convenience [119] with $\lambda'=210\textrm{ nm}$$\lambda'$reference length scale for normalization being the length scale and $\tau' =\lambda'^4 k_BT/ \gamma_s\Omega D_{s} \delta_s$$\tau'$reference time scale for normalization being the time scale. The physical values of the material properties and simulation conditions are listed in Table <ref>.$\rho_{\textrm{mat}}$electrical resistivity of the matrix

[SnAgCu material parameters and experimental conditions adopted in the model]The following SnAgCu material parameters [177] and experimental conditions are adopted in the model.

Atomic volume     $\Omega=2 \times 10^{-29}$ m$^3$

Surface diffusivity     $D_s = 3.99 \times 10^{-10}$ m$^2$/s

Thickness of surface layer     $\delta_s = 2.86 \times 10^{-10}$ m

Surface energy     $\gamma_s = 1.4$ J/m$^2$

Temperature     $T$ = 419 K

Electron charge     $e$ = $1.6 \times 10^{-19}$ c

Boltzmann constant     $k_B$ = $1.38 \times 10^{-23}$ J/K

Electrical resistivity     $\rho_{\textrm{mat}}$ = $1.32 \times 10^{-7}$ $\Omega$m

An explicit finite-difference scheme is utilized for the implementation of the model. The spatial derivatives of the coupled Eqs. (<ref>) and (<ref>) are discretized utilizing a combination of forward, backward, and central differences on a staggered grid, leading to second-order accuracy.
A first-order explicit Euler technique is employed to discretize the temporal derivative. Isolate boundary conditions are imposed on the order parameter $c$, for all directions. The Dirichlet boundary conditions are set at the right and left ends, for the electric potential $\phi$ in the form of prescribed constant value, and are isolate at the top and bottom. The Laplace equation is solved iteratively employing the conjugate gradient scheme.
The termination criterion is selected such that the maximum permissible difference in the electric potential is less than $1 \times 10^{-6}$ at a given position between the current and the previous iteration.

The model is implemented into the in-house software package PACE3DPACE3DParallel Algorithms for Crystal Evolution in 3D version 2.2.0 [178].
In this work, the study is intended to track the complete temporal morphological transition of a circular inclusion.
This may require a extremely large simulation domain of the order of $5000-6000$ grid cells in the direction of inclusion propagation, which increases the computation efforts dramatically.
In the present study, any appreciable change in the field variables $c$ and $\phi$ takes place only in the vicinity of the inclusion, due to the surface diffusion.
The computational cost is restrained within reasonable limits by utilizing the so-called moving window approach [179].
Hence, a simulation starts with a fixed domain size.
When any inclusion part reaches the center of the simulation domain, for each unit cell ($\Delta x$) migration of the inclusion, a grid point is removed from the rear end and is correspondingly added in front of the inclusion.

§ MODEL APPLICABILITY AND LIMITATIONS

Some specific features of the developed model for a transgranular inclusion propagation are discussed as follows:

* The presence of an electric field in the conductors can induce motion of the diffusing metal species in two ways: (1) the direct electrostatic force and (2) the electron-wind force. Firstly, the direct electrostatic force which drives the species towards the negative terminal (cathode). Second, the negatively charged electrons accelerated in the direction of the positive terminal (anode), which are colliding with the species, thus transferring momentum results in the advancement of the species in the direction of electron flow. This contribution is termed as electron wind.  At the operating current densities of interconnects, that is of the order of $10^8-10^{10}$ A/m$^2$, electron-wind force is the dominant of the two. Therefore, within the framework of irreversible thermodynamics, the cross-effect arising due to the interaction between the conducting electrons and the diffusing species is considered. This term has been added phenomenologically in the diffusion equation (<ref>).

* Although the inclusion nucleation, which is due to electromigration, is directly related to the reliability of the interconnects [180], the process of inclusion nucleation is neglected.
Instead, the studies are conducted to specific characteristics associated with the shape evolution of the inclusion.

* The metallic conductor is assumed to be a single crystal, such that the inclusion-grain boundary interactions are neglected.

* The passage of an electric current, through a conductor, produces local heating near sharp corners and bends [181].
Therefore, thermomigration may accompany electromigration. In the presented work, such effects of Joule heating are neglected.
The inclusion shape changes are explored as a result of the competition between the relative magnitude of the electromigration force and the surface capillary force.

* The surface diffusivity at the inclusion-matrix interface is very high in the operating temperature conditions $(T < 500 \textrm{ K})$ [182].
Hence, the surface diffusion is assumed to be the only atomic transport mechanism, and the bulk diffusion is neglected.

* In general, the conductor materials might contain more than a single inclusion [183].
The local electric fields of different inclusions may interact with each other and may alter the overall dynamics of the inclusion evolution.
In the presented dissertation, a study of morphological dynamics of an isolated inclusion is performed neglecting the effect of neighboring inclusions.

§ CONCLUSION

In this chapter, the phase-field method for morphological evolution of inclusion under the external electric field is described. The double-obstacle type of free energy density is considered to represent the temporal evolution of inclusions. The Laplace equation is employed to facilitate the distribution of the electric field in the simulation domain. Besides, numerical strategies considered for model implementation are discussed. The special provision in the mobility equation (<ref>) of the model allows simulating isotropic as well as anisotropic mobilities of inclusions. The results obtained for isotropic inclusions are presented in Chapter <ref>, which also compares sharp-interface theories with numerical results. Thereafter, results on anisotropic inclusions consisting of two-fold symmetry are described in Chapter <ref>, followed by fourfold and sixfold symmetries in Chapter <ref>.

tocchapterIII Results and Discussion:
Phase separation in lithium-ion batteries

Part III

Results and Discussion:
Phase separation in lithium-ion batteries

CHAPTER: SURFACE IRREGULARITIES OF A CATHODE PARTICLE

§ INTRODUCTION

Phase-field literature, which models the Li transport within a single particle [57, 56, 184], is ample. Previous models simulate particle geometry with various simplifications to perform the study in 1D [51, 60] or 2D [185, 160]. There are few exceptional works simulated in 3D, which are limited to regularly shaped geometries, such as spherical particles [58, 186]. The experimental results [187, 188, 189, 190] reveal that the electrode consists of particles that show a shape polydispersity. The effect of particle geometry should be considered to obtain simulated results that are closer to the real electrode particle [191]. Stein and Xu [192] described the influence of particle geometry on the Li concentration profile in ellipsoidal particles. Guo et al. [191] reported that the difference in maximum and minimum Li concentration is larger in elliptical particles than in spherical particles. Chakraborty et al. [193] found that there are important differences between the plastic stretches in cylindrical particles and those in spherical particles, investigated earlier by Cui et al. [194].

The presented work has two objectives. The first objective is to study the insertion dynamics of an irregularly shaped particle, to imitate a real particle in physical existence. The second objective is to consider the effect of various parameters influencing the charge dynamics. The foundation of the model was developed in chapter <ref>. Based on that model, next section <ref> describes employed boundary conditions to study the surface irregularities presented in this chapter. In addition, the theory is validated with commercial multi-physical software COMSOL Multiphysics by considering a particular case from a benchmark study [51]. Thereafter, the simulation results are presented in section <ref> with plausible justification in section <ref>. Afterward, the numerical study is extended in section <ref> to investigate the effect of various material parameters on the transportation mechanism. In section <ref>, the chapter is concluded with some remarks on the results. Some parts of this chapter are published in the journal  Modelling and Simulation in Materials Science and Engineering [59].

§ CONSTANT NEUMANN FLUX CONDITION AT THE PARTICLE SURFACE

For a simplicity, a temporally independent lithium flux $J_n$ is considered at the particle surface of the form,

\begin{equation}
J_n= \boldsymbol{J} \cdot \boldsymbol{n} =     \left\{ \begin{array}{ccl}
\frac{-C R_0}{N_{\textrm{ex}}d \times 3600 } & \mbox{  for}& c_s < 1, \\
0  & \mbox{  for} & c_s = 1.
\end{array}\right. \label{eq:LIBGalvanostaticFlux}
\end{equation}

The parameter $C$$C$C-rate refers to the C-rate, and  x C-rate measures the Li insertion or extraction rate at which complete charging or discharging takes place within 1/x hours [60]. Furthermore, $R_0$$R_0$radius of the reference particle is the radius of the reference particle, $c_s$$c_s$concentration at the surface of the particle is the concentration at the surface of the particle, $d$ defines the dimension of the simulation study, e.g., $d=3$ for a 3D simulation study and the characteristic extension$N_{\textrm{ex}}$characteristic extension

\begin{equation}\label{Eq:LIBGalvanostaticCurvature}
N_{\textrm{ex}} =     \left\{ \begin{array}{rcl}
\frac{\mathrm{Particle\, circumference \,} \times \mathrm{ \,Reference\, particle\, area}}{\mathrm{Particle \,area \,} \times \mathrm{\, Reference\, particle\, circumference}}& \mbox{  for the}& \mathrm{2D\, case,} \\
\frac{\mathrm{Particle\, surface\, area \,} \times \mathrm{ \,Reference\, particle\, volume}}{\mathrm{Particle\, volume \,} \times \mathrm{\, Reference\, particle\, surface\, area}}  & \mbox{  for the} & \mathrm{3D \,case.}
\end{array}\right.
\end{equation}

The constant Neumann flux boundary condition in Eq. (<ref>) on the arbitrary geometrical shapes is implemented through the evolution equation (<ref>) for the Li diffusion. In addition, homogeneous nucleation from the surface is considered in the form of Eq. (<ref>).  For the 2D (or 3D) simulation study, a circular (or spherical) particle with radius R$_0$ = 1.0 $\mu$m is considered as a reference particle along with other parameter set is provided in Table <ref>. A study is conducted to observe the Li diffusion from the surface of the particles, which are reported here.

§.§ Simulation domain setup

The general boundary condition of arbitrary geometrical shapes is implemented through the smoothed boundary method [160, 161], in which a domain parameter, $\psi_1$, is registered to differentiate and interpolate between different phases. The simulation is initialized with the domain parameter $\psi_1$, continuously defined in the entire domain. The field quantity $\psi_1$ takes the value 1.0 for the bulk cathodic particle, 0.0 for the electrolyte and varies smoothly in the boundary between the particle and the electrolyte. The specific features of a diffuse interface between the particle and the electrolyte are responsible for the configuration of the surface normals, as shown in Figure <ref>. As a consequence, the particle surface normal is exploited into the flux boundary condition given in Eq. (<ref>).

In the present chapter, particle morphology is considered to be constant during lithiation. Therefore, the simulations are performed in two independent steps. In the first step, the particle is formed from the Allen-Cahn equation presented in Appendix <ref>. The resultant domain encompasses a smooth interface between the cathode particle and the electrolyte. In the second step, the smooth surface is exploited to implicitly incorporate the flux boundary condition in the Cahn-Hilliard equation (<ref>), in which lithiation inside the particle obtained from the first step is considered.

The model simulation of the species flux infiltrates from the lithium abundant electrolyte to the cathode particle, through the exterior surface. A simulator is developed to obtain the numerical solution of the compound partial differential equation of the Cahn-Hilliard equation combined with a smoothed boundary method. By employing standard message passing interface (MPI)MPImessage passing interface concepts, a parallel three-dimensional solver for large-scale computations on a rectangular mesh is realized as described in section <ref>.

It is important to note that owing to the smoothed boundary method, the electrolyte is considered as LMO material as well. Even though the evolution equation (<ref>) is solved everywhere in the domain, i.e. in the region of the particle ($\psi_1=1$) as well as in the electrolyte ($\psi_1=0$), the specific manipulation in the smoothed boundary method motivates no concentration accumulation at $\psi_1=0.0$ (i.e., electrolyte). As a consequence, the influence of the electrolyte material parameter on the concentration accumulation in the electrolyte is negligible.

[Material properties of LiMn$_2$O$_4$ cathode and operational conditions]Material properties of LiMn$_2$O$_4$ cathodic intercalation materials [51] and operational conditions.

Parameter     Symbol     Value      Unit

Diffusion coefficient     D$_1$     7.08 $\times$ 10$^{-15}$     m$^2$/s

Length scale     $L$     1     $\mu$m

Gradient energy coefficient     $\kappa$     7 $\times$ 10$^{-18}$     m$^2$

Regular solution parameter     $\alpha_1'$     -0.1     -

$\alpha_1''$     2.6     -

Temperature     $T_{\mathrm{ref}}$     300     K

§.§ Validation of the presented approach

[Effect of the dimensionless cell width and the interface thickness coefficient on the SOC value]Effect of the dimensionless cell width and the interface thickness coefficient on the SOC value: Comparison of the present approach with results obtained by the finite element software COMSOL Multiphysics®.

The simulation studies are performed for LMOLMOlithium manganese oxide, LiMn$_2$O$_4$ cathodic particle at a lithiation of default 1 C-rate discharging condition. Therefore, the inflow of the lithium flux fills up the particle. The natural logarithmic terms $\mathrm{ln}({c})$ (or $\mathrm{ln}(1-{c})$) in the chemical-free energy formulation are singular at a concentration of ${c}=$ 0.0 (or 1.0). To avoid the uncertainty, the particles are filled with a concentration of a very small positive fraction, e.g., ${c_i}=$ 0.03, unless stated otherwise. Similarly, the simulations stop when any of the positions reaches ${c}=$ 0.999, instead of a fully lithiated state, inside the particle.

Yu et al. [161] presented a smoothed boundary method based on a finite difference scheme for various applications, such as the mechanical equilibrium equation, phase transformation with the presence of additional boundaries and the diffusion equation with Neumann and Dirichlet boundary conditions. In this section, the presented approach in section <ref> is validated by considering a particular case in a phase-separating particle. During the insertion of the lithium species, a concentration gradient develops inside the particle from the higher concentration at the surface to the lower concentration at the center [195]. As a consequence, the surface reaches the threshold value of $c=0.999$ to stop the simulation before the bulk region of the particle. In turn, the simulations stop well before the particle attains a SOCstate of chargestate of charge (SOC = $100/{V_\Omega} \int_{V_\Omega} \psi {c})$ of 100% for 1 C-rate species inflow. Therefore, the latest SOC attained by different values of simulation cell width of spherical particles are compared, with the 1D computations performed in Ref. [135] by using the software COMSOL Multiphysics®, which is based on the finite element method. The results obtained from the presented model are in agreement with the study of Huttin and Kamlah [51]. The comparison suggests that at most a dimensionless cell width of $\Delta \boldsymbol{x} = 5 \times 10^{-2}$ is acceptable, to obtain a sufficient convergence of the simulation results, see Figure <ref>. Moreover, the calculated concentration value deviates further from the benchmark value as the interface thickness increases. The results ensure convergence for a thin interface thickness, which is in agreement with literature [161]. Therefore, the simulation study is carried out at $\Delta \boldsymbol{x} = 1.67 \times 10^{-2}$ for 3D3Dthree-dimensional and at $1 \times 10^{-2}$ for 2D2Dtwo-dimensional studies with the interface thickness coefficient$\epsilon$interface thickness coefficient ${\epsilon_p}/R_0=0.05$, which are confined to the convergence of the smoothed boundary method.

§ CONCENTRATION PROFILES

[Temporal evolution of isolated 2D elliptical particles with an initial concentration of ${c_i}=$ 0.03.]Temporal evolution corresponding to a 10%, 60%, and 90% state of charge, including an intermediate concentration profile at the initiation of a phase separation inside 2D elliptical particles, with aspect ratios $\gamma=$ 1, 2, 3, 4, and 5. The gold and magenta colors respectively correlate to the Li-rich and Li-poor phases. The particles are filled with an initial concentration of ${c_i}=$ 0.03.

[Concentration profiles of isolated 2D elliptical particles with an initial concentration of ${c_i}=$ 0.26.]Concentration profiles corresponding to 30%, 50%, 70%, and 90% state of charge inside 2D elliptical particles, with aspect ratios $\gamma=$ 1, 2, 3, 4, and 5. The gold and magenta colors correlate to the Li-rich and Li-poor phases, respectively. The particles are filled with an initial concentration of ${c_i}=$ 0.26.

To begin with, a 2D simulation study with a varying aspect ratio of ellipses of equal area is considered. The ellipses with the aspect ratios $\gamma$$\gamma$aspect ratio of ellipse = 1, 2, 3, 4, and 5 with major semiaxes are 1.00, 1.42, 1.74, 2.00, and 2.25 $\mu$m, respectively in Figure <ref> with initial filling ${c_i} =0.03$$c_i$initial concentration filling and in Figure <ref> with ${c_i} =0.26$. The first column in Figure <ref> shows that a 10% state of charge (SOC) corresponds to a homogeneous concentration profile. It is evident that the phase separation initiates in the particles with a higher aspect ratio, prior to the lower ones. For instance, phase segregation starts at 20.38% SOC for $\gamma=1$, while at 18.56% SOC for $\gamma=5$. Furthermore, the circular particle (i.e., $\gamma=1$) evolves with a continuous phase separation across the interface, while two Li-rich islands form at the maximum curvature points in the ellipses (i.e., $\gamma=2, 3, 4,$ and 5). The Li-rich islands continue to grow further at the particle surface. In Figure <ref>, beyond 60% SOC, the Li-poor phase is surrounded by the Li-rich phase in all particles, meaning they eventually follow the “core-shell model", and the core shrinks as the shell enlarges with the evolution of time, which is called the “shrinking-core model". Contrary to ${c_i} =0.03$ cases, there are no distinct Li-rich islands are observed for ${c_i} =0.26$, instead a continuous region of Li-rich phase completely surrounds the Li-poor phase for each particles as shown in Figure <ref>.

[Concentration profiles of 2D irregularly shaped particle] Concentration profiles of an irregularly shaped particle, when the SOC ranges from approximately 19% to 90%. The images correspond to 19%, 36%, 77% and 90%. The red color, $\psi_1$ = 1.0, corresponds to the cathodic particles, and the blue color, $\psi_1$ = 0.0, corresponds to the electrolyte in (a). The gold color represents the lithium-rich phase, and the magenta color inside the lithium-rich phase visualizes the lithium-poor phase from (b) onwards.

Figure <ref> shows the evolution history of a 2D irregularly shaped particle. It is evident that the Li-rich phase initiates at the maximum curvature regions, as indicated by the three gold color islands in Figure <ref>(b). These islands grow and, as a consequence, merge to form a single continent. Simultaneously, the Li-poor phase takes a curvaceous shape, followed by an ellipse and eventually transforms into a circle as it shrinks, which can be seen in Figure <ref>(b), (c) and (d) respectively.

A 2D simulation study of five elliptical particles simultaneously immersed in the electrolyte, see Figure <ref>. The particles are filled with an initial concentration of ${c_i}=$ 0.26. The particle aspect ratio varies from $\gamma = $ 1, 2, 3, 4 and 5 with 1.0, 1.54, 2.13, 2.72, and 3.35 $\mu$m major semiaxes and 1.00, 0.77, 0.71, 0.68, and 0.67 $\mu$m minor semiaxes ellipses, respectively. The circumference to area ratio of the particles is fixed. Therefore, $N_{\textrm{ex}}$ = 1 for all particles according to Eq. (<ref>) and consequently, the particles are subject to identical applied flux and the particles evolve with identical SOC. The evolution pattern follows the widely known “core-shell model."

[Concentration profiles of five elliptical particles simultaneously immersed in an electrolyte with initial concentration $c_i = 0.26$]Concentration profiles
of five particles, simultaneously immersed in an electrolyte, when the SOC ranges from approximately 30% to 90%, displayed with an increment of 20%, until it reaches 96.37% in the last frame (f). In the entire domain, the domain parameter $\psi_1$ is continuously defined. The yellow color, $\psi_1$ = 1.0, corresponds to the cathodic particles and the blue color, $\psi_1$ = 0.0, corresponds to the electrolyte in (a). The gold color represents the lithium-rich phase, and the magenta color inside the lithium-rich phase visualizes to the lithium-poor phase from (b) onwards.

Figure <ref> shows the evolution history of a 3D particle with an irregular shape. The Li-rich islands, which are marked in red color, develop at the convex surface, which corresponds to the maximum curvature of the particle. The islands progress to form a single continent at 40% SOC and eventually, the core shrinks with time.

For various initial concentrations, a comprehensive study is performed as a function of aspect ratio $\gamma$ for a SOC at which phase separation is expected, SOC$_\mathrm{PS}$state of charge at the onset of phase separation$\mathrm{SOC}_\mathrm{PS}$, see Figure <ref>(a). The curve fit suggests the empirical relation,

\begin{equation}\label{eq:SOCvsGamma}
\mathrm{SOC}_\mathrm{PS} = (-6.8  {c}_i^2 +20) \gamma^{(1.3  {c}_i^2-0.061)} \textrm{ for } 0.00<{c}_i<0.20,
\end{equation}

where, SOC$_\mathrm{PS}$ is the SOC at the onset  of phase separation. It is evident that increasing the initial concentration reduces the difference in the values of SOC$_\mathrm{PS}$ between the higher and the lower aspect ratio particles continuously. As a consequence, at a certain initial concentration, namely ${c_i} = 0.26$, phase separation starts in all particles simultaneously as seen in Figure <ref>.

[Effect of initial concentration ${c}_i$ on the phase separation dynamics]The SOC (or duration) at which phase separation initiates, SOC$_\mathrm{PS}$ in (a) (or $t_{\mathrm{PS}}$$t_\mathrm{PS}$time at the onset of phase separation in (c)) and completesSOC$_\mathrm{PE}$state of charge attained before the end of phase separation SOC$_\mathrm{PE}$ in (b) (or $t_{\mathrm{PE}}$$t_\mathrm{PE}$time attained at the end of phase separation in (d)) are plotted as a function of the aspect ratio of the 2D elliptical particles for the initial concentrations ${c_i} =$ 0.03, 0.05, 0.11 and 0.13. The dotted lines in (a) are plots of Eq. (<ref>), while in (b), (c) and (d) are line fits for respective values of ${c}_i$.

[Concentration profiles of 3D irregularly shaped particle]Concentration profiles of an irregularly shaped particle, when the SOC ranges from approximately 20% to 90%. The images correspond to 20%, 40%, 60%, 80% and 90%. The opaque silver colored figure in (a) indicates the contour of the particle surface at the value $\psi_1$ = 0.5 of the domain parameter, and the particle from the upper section is clipped for the sake of a better visual sight from (b) onwards. The red color is a contour of ${c}=0.80$ and it represents the lithium-rich phase. The blue color is a contour of ${c} = 0.22$ and the blue color cavity inside the lithium-rich phase corresponds to the lithium-poor phase.

The basic features of the charge dynamics and its dependence on the various model parameters are explained with feasible reasons in the following section.

[Phase diagram indicating various zones characterized by the miscibility and the spinodal gaps]The phase diagram shows various zones, which are characterized by the miscibility gap and the spinodal gap. Zones I, II and III respectively correspond to homogeneous states, nucleation states, and phase-separated states. The solid blue line in the bottom graph is a chemical-free energy density, $f_1^{\mathrm{chem}}$ at absolute temperature, $T=300$ K. The projections of the miscibility gap, on the free energy density curve, indicate the equilibrium conditions expressed in Eqs. (<ref>) and (<ref>). The projections of the spinodal gap, on the free energy density curve, indicate zero curvature points. The dotted red line in the bottom graph corresponds to the minimum energy path, widely known as Maxwell construction.

§ MISCIBILITY GAP

The local concentration characterizes a condition of the particles, which corresponds to the zones I, II and III in Figure <ref>, which determine a phase-segregated state or a homogeneous state during intercalation. In the phase equilibrium diagram of the LMO particle, these zones are related to spinodal and miscibility gaps. In Figure <ref>, the spinodal gap is the locus of points where the curvature of the chemical free energy density changes its sign, i.e. from convex to concave, at point B and from concave to convex, at point C. Let $c^L$$c^L$lower concentration level of a phase-separated state and $c^H$$c^H$higher concentration level of a phase-separated state be the two concentration levels of a phase-segregated state for the Li-rich and the Li-poor phases, respectively. In electrode particle ($\psi_1=1$), the miscibility gap is obtained by respectively equating the chemical potential $\mu_1^{\mathrm{chem}}$ and the grand-chemical potential $G_a^{\mathrm{chem}}$grand-chemical potential of phase $a$$G_1^{\mathrm{chem}}=f_1^{\mathrm{chem}} - c \mu_1^{\mathrm{chem}}$ for both concentrations, $c^L$ and $c^H$, as

\begin{equation}
\mu_1^{\mathrm{chem}}(c^L) = \mu_1^{\mathrm{chem}}(c^H), \label{Eq:potential}
\end{equation}

\begin{equation}
f_1^{\mathrm{chem}}(c^L) - c^L \mu_1^{\mathrm{chem}}(c^L) =f_1^{\mathrm{chem}}(c^H) - c^H \mu_1^{\mathrm{chem}}(c^H).\label{Eq:grandPotential}
\end{equation}

The simulation study is carried out at an absolute temperature of $T=$ 300 K. The corresponding chemical free energy density plot in Figure <ref> represents the double-well structure, which entails the presence of a concavity between the two wells, and hence the presence of miscibility and spinodal gaps. The concentration of the miscibility gap fulfills the equilibrium conditions expressed in Eqs. (<ref>) and (<ref>), which are marked by the points A and D. Between these points, the free energy can be decreased by decomposing the homogeneous state of higher energy to a phase-segregated state of lower energy. The spinodal gap concentrations between the zero curvature points of the free energy curve, marked as points B and C, render the homogeneous distribution states unstable and force the particle into phase segregation. In Figure <ref>, the red dotted line refers to the common tangent construction and represents the minimum energy line.

With the application of flux, a concentration gradient is established from the surface to the center of the particle. Therefore, the phase segregation initiates from the surface of the particle. The generation and increase in the pre-existing phase boundary between the Li-rich and the Li-poor phases are energetically costly, and the system follows a path with the minimum phase boundary possible. The increase in the phase boundary portion leads to an energy penalty. Thus, the equilibrium states tend towards a minimum phase boundary inside the particle, following the evolution in the direction of the minimum free energy path. As a consequence, the inner phase shrinks by assuming an irregular shape with a curved boundary, at the beginning, which then transforms into an ellipse, and eventually into a circle, as shown in Figure <ref>.

In Figure <ref>, the uniform curvature along the surface of the circular particle ($\gamma$ = 1) results in a continuous phase separation across the particle surface. In the case of $\gamma$ = 2, 3, 4 and 5, two Li-rich islands develop at the maximum curvature points. The likely reason behind the formation of the islands is that the local concentration at these positions reaches the spinodal zone first and causes a local phase separation. In an ellipse, the largest and shortest distances from the center of the ellipse to the surface correspond to the major and minor semiaxes. There, the maximum and minimum curvature points are also located. In the vicinity of the maximum curvature points, the particles have a larger circumference to the adjacent particle area, i.e., a larger local circumference to area ratio, which is in contrast to the minimum curvature points. As a consequence, a higher concentration builds up at the surroundings of the endpoints of the major semiaxis, in contrast to the minor semiaxis, which is in agreement with the literature [191, 192, 196]. This comparison suggests that more sites are available for the applied flux, compared to the sites that host lithium species near the major semiaxis. This triggers a phase separation process during the insertion as seen in Figure <ref>.

[Maximum and minimum concentration, along the surface of the particles, for the initial concentrations ${c}_i = 0.03$ and $0.13$ at 18% SOC]The maximum and minimum concentration, along the surface of the particles, for the initial concentrations ${c}_i = 0.03$ and $0.13$, at 18% SOC. The maximum concentration observed in a region with maximum curvature, while the minimum concentration is observed in a region with minimum curvature in ellipses.

Aspect ratio, $\gamma$     2 c Maximum concentration, ${c}$     2 c Minimum concentration, ${c}$

${c_i}=0.03$     ${c_i}=0.13$     ${c_i}=0.03$    ${c_i}=0.13$

1     0.200     0.201     0.200     0.201

2     0.216     0.214     0.187     0.188

3     0.224     0.220     0.180     0.182

4     0.231     0.223     0.175     0.179

5     0.239     0.225     0.172     0.177

The difference between the maximum and minimum concentration across the surface of the particle increases with the aspect ratio, as can be seen in Table <ref>. Therefore, the maximum concentration point reaches the spinodal gap with less SOC for particles with a higher aspect ratio. Hence, the SOC required to initiate a phase separation decreases with the aspect ratio. Similar to Figure <ref>, the higher curvature region accumulates more lithium, and arrives at the spinodal region, prior to the lower curvature region. Therefore, the phase separation starts at the corners and edges which correspond to the higher curvature region. In the meantime, it can be observed that the particle surfaces close to the minor semiaxis accumulate further species. The Li-rich islands continue to grow and eventually, the islands merge to form a single continent.

§ INFLUENCE OF VARIOUS PARAMETERS ON PHASE SEPARATION MECHANISM

§.§ Effect of initial concentration

Figure <ref>(a) shows the increment in SOC$_{\mathrm{PS}}$ for higher aspect ratio particles with a rise in initial concentration ${c}_i$. An explanation is the decrease in the difference between a maximum and a minimum concentration at higher ${c}_i$, see Table <ref>. Furthermore, an increase in ${c}_i$ results in less concentration required locally to reach the spinodal zone, and drive phase separation. If the initial concentration ${c}_i$, in all particles of Figure <ref> falls in the spinodal zone, phase separation is initiated everywhere at the surface of the particles simultaneously, regardless of the curvature variation across the particle surface. Furthermore, the particle with the higher ${c_i}$ requires less time to reach the same SOC compared to the lower. Therefore, the duration required to observe the onset of phase separation (i.e., $t_{\mathrm{PS}}$) decreases for higher initial concentration, as reported in Figure <ref>(c).

Although the initiation of the phase separation process is affected by the initial concentration of ${c_i}$, the completion of phase separation is independent, as depicted in Figure <ref>(b). During the spinodal decomposition, the particle attains two concentration levels (i.e., Li-rich and Li-poor phases), irrespective of the initial concentrations. In other words, the phase separation process hinders the effect of initial concentration. Therefore, the effect of different initial concentrations on the SOC required to complete the phase separation ( i.e., SOC$_\mathrm{PE}$) is not observed. In addition, the effect of variation of ${c_i}$ on the duration until the phase separation observed $t_{\mathrm{PE}}$ is reported in Figure <ref>(d). Obviously, the increase in ${c_i}$ shows additional leverage to SOC, which does not require to accumulate during operation. As a consequence, the increase in ${c_i}$ leads to a decrease in the time to complete the phase separation $t_{\mathrm{PE}}$.

[Variation of C-rate and mobility as a function of the aspect ratio of the elliptical particles]Variation of C-rate as a function of the aspect ratio of the elliptical particles for SOC at which the onset of phase separation starts  SOC$_{\mathrm{PS}}$  in (a) and ends SOC$_{\mathrm{PE}}$ in (b). The variation of mobility $\hat{M}$ as a function of SOC$_{\mathrm{PS}}$ in (c) and SOC$_{\mathrm{PE}}$ in (d).

§.§ Effect of C-rate and mobility

Ideally in the quasi-static condition (i.e., C-rate $C\ll 1$) [51], the phase separation starts when the average concentration reaches the spinodal point B and ends outside the miscibility gap  (point D in Figure <ref>). In fact, the local concentrations at all positions are almost equal to the average concentration in quasi-static conditions. However, higher values of C-rate stimulate a steeper concentration gradient from the surface to the center. Therefore, a local concentration at the surface reaches the spinodal point before the average concentration. The difference between these two increases with flux. As a result, SOC required to observe the onset of phase separation SOC$_{\mathrm{PS}}$ decreases for higher flux rates, which is depicted in Figure <ref>(a). Contrarily, the phase separation ends when the Li-poor phase is completely depleted by the Li-rich phase. The slope of the concentration gradient increases with the C-rate. As a consequence, the presence of gradient penalizes the SOC required to end the phase separation SOC$_{\mathrm{PE}}$ on a higher side, which can be seen in Figure <ref>(b).

On the one hand, mobility can be related to the inherent ability of the system to homogenize the concentration distribution by transferring species in order to eliminate any gradients developed during operation. On the other hand, C-rate enhances the concentration gradient at the surface by depositing the flux of lithium species. The interplay between the deposition and transfer rates characterizes the equilibrium states. At the elevated mobility, the species transfer rate is more prominent compared to the species deposition. As a consequence, the species disperse immediately after the deposition at the surface. Thus, a particle requires higher SOC to observe phase separation at higher mobility as depicted in Figure <ref>(c). Furthermore, the slope of concentration decreases at the agile mobility, which results in decreased SOC$_{\mathrm{PE}}$ (see Figure <ref>(d)). Conversely, the lower mobility provides comparatively enough time for the accumulation of the species, which are deposited at the particle surface. Hence, SOC$_{\mathrm{PS}}$ decreases and SOC$_{\mathrm{PE}}$ increases for sluggish mobility as an effect from the increased slope of the concentration gradient.

A very high mobility (or very low C-rate) corresponds to the steady-state, in which phase separation starts at spinodal point B in Figure <ref>, for all particles irrespective of curvature. In other words, the curvature effects are suppressed in steady-state. As the mobility decreases (or C-rate increases), the curvature effects become prominent as a result of the development of concentration gradients. For instance, the increase in the difference between SOC$_{\mathrm{PS}}$ for the $1.0\hat{M}$ and $1.2\hat{M}$ with aspect ratio $\gamma$ is plausible due to the particle curvature.

[Effect of temperature $T/T_{\mathrm{ref}}$ and free energy parameter $\alpha_1''$ on free energy density plots]The effect of variation in (a) temperature $T/T_{\mathrm{ref}}$ and (b) free energy parameter $\alpha_1''$ on free energy density plots. The points A and D represent free energy local minima obtained from Eqs. (<ref>) and (<ref>), while B and C are spinodal points where the curvature of the plot changes its sign.

§.§ Effect of temperature and free energy parameter $\alpha_1''$

[Variation of the temperature and free energy parameter as a function of the aspect ratio of the elliptical particles]Variation of the temperature $T/T_{\mathrm{ref}}$ as a function of the aspect ratio of the elliptical particles for SOC at which the onset of phase separation starts  SOC$_{\mathrm{PS}}$ in (a) and ends SOC$_{\mathrm{PE}}$ in (b). The variation of free energy parameter $\alpha_1''$ as a function of SOC$_{\mathrm{PS}}$ in (c) and SOC$_{\mathrm{PE}}$ in (d).

During operating conditions, the battery is seldom subjected to a constant temperature. Thus, understanding the effect of variation in operating temperature on the phase separation dynamics has technological implications. The free energy density for three different values of the temperature $T/T_{\mathrm{ref}}$ is plotted in Figure <ref>(a). The particle, which is initially filled with concentration value of ${c_i}=0.03$, starts accumulating species during insertion and consequently increases in local concentration. Note that, if concentration at a particular position of the particle reaches spinodal zone (i.e., point B), it decomposes into two concentration levels (i.e., lower concentration point A and higher concentration point D). The rise in temperature shifts the spinodal point B on the higher side (see Figure <ref>(a)). In other words, the system requires a higher value of concentration so that any of the positions lie in the spinodal zone at the elevated temperature. Therefore, SOC required to observe the phase separation SOC$_\mathrm{PS}$ increases with temperature $T / T_{\mathrm{ref}}$, which can be seen in Figure <ref>(a). On the contrary, the rise in temperature shortens the miscibility gap, i.e., a shift in point D towards the lower concentration values (see Figure <ref>(a)). As the point D is closely related to SOC$_\mathrm{PE}$, the effect of a variation of the temperature $T/T_{\mathrm{ref}}$ is justified in Figure <ref>(b).

Figures <ref>(c) and (d) show the effect of variation of the free energy parameter $\alpha_1''$. In Figure <ref>(b), it can be seen that the miscibility and spinodal gaps enlarge with an increase in $\alpha_1''$, the behavior is opposite to the temperature $T/T_{\mathrm{ref}}$.  Therefore, SOC$_\mathrm{PS}$ decreases and SOC$_\mathrm{PE}$ increases for higher $\alpha_1''$, which shows an opposite tendency compared to the temperature $T/T_{\mathrm{ref}}$.

Finally, Figure <ref> and Figure <ref> exhibits the capability of the method to simulate multiple particles simultaneously and three-dimensional particles with an irregular shape immersed in an electrolyte. As the particle boundary condition is implicitly defined in the smoothed boundary method, this method can be applied to particles with almost any geometry. Hence, this technique is very powerful and convenient to solve differential equations in complex geometries with complex boundary conditions that are often difficult to mesh.

§ CONCLUSION

A phase-field study is performed to simulate phase separation during the insertion process. The presented model employs two coexisting phases, which is validated with a benchmark. In addition, the increase in mesh resolution conveys the convergence of the presented method. In the finite-difference framework, the smoothed boundary method enables the modeling of particles with almost any geometry. The effects of the particle geometry on the concentration evolution are explored numerically. For a constant flux boundary condition at the particle surface, the free energy density and the chemical potential are discussed in detail. Based on that, the spinodal and miscibility gaps are estimated.

The simulation results show that the growth of the phase boundary between the Li-poor and Li-rich phases, inside the particle, leads to an energy penalty. As a consequence, the morphological evolution of the concentration profile suggested a minimum phase boundary pathway. With the application of flux, the particle surface with a higher curvature preferentially accumulates lithium species. Therefore, phase segregation starts in the vicinity of the regions with a higher curvature. The reason for the inhomogeneous phase separation across the particle surface is explained in detail, by means of miscibility and spinodal gaps. Furthermore, the elliptical particle with a higher aspect ratio is subjected to the onset of the phase separation, prior to the lower ones.

Finally, it is evident from Figure <ref> that the smoothed boundary method can be applicable to multiple particles simultaneously. However, the significance of the simulation involving multiple particles is restricted due to the boundary conditions are applied on the particle surfaces. Therefore, the influence of neighboring particles is limited. In fact, the transportation of species in any particle is independent of the other particles and therefore could be simulated separately of the others, as shown in Figure <ref>. Hence, the concentration profiles are identical for the densely packed particles and the sparsely packed, which should not be the case generally. This provokes an extension of the simulation study to identify the significance of the microstructural properties of the electrodes consists of many particles. Therefore, these issues are addressed in the next chapter.

CHAPTER: MORPHOLOGICAL DESCRIPTORS IN MULTIPLE PARTICLE POROUS ELECTRODES

High-performance batteries necessitate electrodes of superior active material and more importantly, optimized porous microstructures [66, 62]. Recent advancements in manufacturing facilities provide great control over the particle size and the complex-tortuous microstructure of electrodes [197]. Thus, electrode engineering allows producing high-density electrodes without compromising its rate capability. The physical microstructure of the electrode has a direct impact on battery performance measures such as cyclic stability, intercalation rates, and power density amongst others [198]. Therefore, it is important to quantify the parameters that control electrode microstructure, such as particle shape, size, porosity, and tortuosity. For instance, an experimental study [199] reported that the rate performance of the graphite electrode to be a function of particle size and porosity of the electrode. In addition, images from X-ray tomography suggest that the microstructural inhomogeneities largely influence the direction of lithium transport [200, 201].

Systematic numerical studies of porous electrodes may complement the scientific understanding of the hierarchical microstructures and help to optimize the process parameters efficiently and economically [202, 61, 203]. Despite some theoretical and numerical studies are performed to account essential features of phase-separating porous electrodes [63, 37, 64, 65], the literature providing exclusive information on the microstructural properties of the electrodes is scarce. For instance, Ferguson and Bazant [204] showed a change in effective diffusivity relates to the system porosity and tortuosity. Besides, Orvananos et al. [67] demonstrated the impact of electrode architecture in terms of ionic and electronic connectivities by considering interactions between active particles. Furthermore, Vasileiadis et al. [44] investigated that the tortuosity, particle sizes, porosity, and electrode thickness influence the capacity of Li$_4$Ti$_5$O$_{12}$ (LTO)LTOlithium titanate, Li$_4$Ti$_5$O$_{12}$ electrodes as a function of C-rate. However, most of the studies utilize an approximation by the Bruggeman relation [205] for the calculation of tortuosity, instead of estimating the tortuosity based on the electrode microstructure under investigation [206, 207, 208]. In addition, these porous electrode phase-field models are focused on LiFePO$_4$ or Li$_4$Ti$_5$O$_{12}$ electrode materials. Relatively few models have been developed for LMO materials [51, 55, 59, 60], which are limited to single particles. Therefore, this chapter aims to provide crucial measures to define the electrode microstructures and the relation of those attributes to the performance of LMO electrodes. Additionally, the microstructural properties, specifically the tortuosity, are estimated numerically to provide a better understanding of the transport mechanism [206].

In the present chapter, the transportation rates of porous electrode structures with the two-phase coexistence in LMO particles is investigated. In this initial effort, neglecting the reaction kinetics at the electrode-electrolyte interface, the focus is based on the diffusional properties of the electrode and the electrolyte. In addition, an explicit treatment of mechanical effects due to misfit strains at the phase boundary is avoided. Instead, the emphasis is placed on the influence of various microstructural properties such as particle sizes, porosity, and tortuosity of electrode consists of ellipsoid-like particles with a focus on the transport mechanism. The chapter organized as follows: Section <ref> illustrates considered simulation setup to investigate the electrode microstructures. Few typical cases of insertion dynamics are discussed in section <ref> with highlighting the underlying mechanism. Also, the obtained results are analyzed and compared with known relations. Thereafter, a numerical study of various morphology descriptors is performed in section <ref>. Finally, the chapter is concluded by a brief discussion on the applicability of the presented results in Section <ref>. Some parts of this chapter are submitted for publication as a journal article.

§ SIMULATION SETUP

The generic model of phase separation for multiphase systems of LMO material is described in chapter <ref>. In the present chapter, the evolution equation (<ref>) is solved numerically with considering a constant concentration boundary condition (<ref>), as shown in Figure <ref>.  As a demonstration of the model, the results are obtained for considering only two phases in the present chapter, the electrode and the electrolyte with highlighting the broad applicability. Particularly, the simulation results are presented for propylene carbonate +1$M$ LiClO$_4$ (lithium perchlorate) as an electrolyte and Li$_x$Mn$_2$O$_4$ as an electrode. The utilized values of parameters for the simulation study are expressed in Table <ref>. The characterization of the electrode microstructure consists of several particles is described in Appendix <ref>.

The simulation setup consists of cathode particles and the electrolyte. The respective regions are filled with their equilibrium concentration values to avoid any self-generated driving force to the diffusing species, which may hinder the effect of external lithium driving force (constant Dirichlet concentration condition). In the present work, the study is intended to focus the insertion of lithium species. Therefore, the electrode region is filled with lower equilibrium value ${c_i}(i, j, k, {t}=0)=0.12$, while the electrolyte region with its respective equilibrium ${c_i}(i, j, k, {t}=0)=0.50$. A higher concentration value is prescribed at the boundary ${c}_{ps}(i={N}_x, j, k,{t})$(= 0.57 in the present study) to provide a driving force to the diffusing species. The simulation domain is observed by monitoring the state of charge (SOC) frequently, once there are no changes observed for a reasonable time then the simulations can be stopped. Note that the same model can be utilized to study the lithium extraction process. Then, the minute changes are need to be done in terms of electrode regions filled with higher equilibrium value and the prescribed constant concentration boundary condition should be lower than 0.5.

[Simulation conditions and material properties of propylene carbonate +1$M$ LiClO$_4$ electrolyte and Li$_x$Mn$_2$O$_4$ electrode.]Simulation conditions and material properties of propylene carbonate +1$M$ LiClO$_4$ electrolyte and Li$_x$Mn$_2$O$_4$ electrode.

Parameter   Symbol   Value   Unit

Reference Temperature   T$_{\textrm{ref}}$   300[59]   K

Absolute Temperature   T   300[59]   K

Reference length scale   $L$   55.5    nm

Electrode Parameters:

Diffusion coefficient   D$_1$   $1\times 10^{-13}$[209]   m$^2$  s$^{-1}$

Gradient energy coefficient   $\kappa_1$   $7\times 10^{-18}$[51]   m$^2$

Regular solution parameters   $\alpha'_{1}$   2.5[51]   -

$\alpha''_{1}$   -5.2[51]   -

Electrolyte Parameters:

Diffusion coefficient   D$_2$   $2.58\times 10^{-10}$[209]   m$^2$  s$^{-1}$

Gradient energy coefficient   $\kappa_2$   $0.0$ (Appendix <ref>)   m$^2$

Regular solution parameter   $\alpha''_{2}$   0.0 (Appendix <ref>)   -

§ BASIC FEATURES OF THE MODEL

§.§ Concentration profiles for porous and planar electrodes

[Representative cases of insertion of species in 2D porous electrode under a constant concentration boundary condition at the separator]Representative cases of insertion of species in a porous electrode with $\unit[11.1]{\mu m}$ depth along with current collector under a constant concentration boundary condition at the separator in two-dimensional space. (a) shows insertion in particles aligned to the species flow, i.e. perpendicular to the current collector, while (b) shows particle parallel to the current collector. The leftmost column demonstrates simulation geometry, while the other four columns depict phase-separated morphologies for various time-steps and SOCs.

Figure <ref> shows insertion in particles aligned perpendicular to the current collector in (a) and parallel to the current collector in (b) for equal porosity and particle size. Firstly, the onset of the development of the Li-rich phase is observed near the separator initially, where the constant concentration boundary condition is employed as evident in (a2) and (b2). On the one hand, the coexistence of two-phases can be observed for a greater depth of electrode in the particles perpendicular to the current collector in (a3). On the other hand, even after a longer period, the coexistence of two-phases is observed for a shorter length span of the electrode in (b3), along with completely intercalated and deintercalated particles near the separator and the current collector respectively. The former case corresponds to more concurrent intercalation, while the latter case rather follows the particle-by-particle intercalation [37].

Within a single particle, due to the diffusing species conveniently cover the depth of almost the complete electrode region in the former case, phase separation initiates nearly from the whole surface of a particle, which is in contrast to the latter where phase separation observed mainly at the regions oriented toward the separator as shown in (b4) compared to (a4). Therefore, the core-shell type phase separation within a single particle is more prominent in the former case as shown in (a4), while a phase-separated traveling front from the separator to the current collector is apparent in the latter case in (b4) for the porous electrode.

[Average of concentration ${c}$ in porous and planar electrodes as a function of the depth of the electrode from the separator]Average of concentration ${c}$ as a function of the depth of the electrode from the separator for various SOC values. The concentration profiles of the porous electrode (Figure <ref>(a) and (b) correspond to (a) and (b) of the present Figure respectively) is compared with the planar single electrode particle model.

The porous structure of an electrode is of interest in battery applications due to high surface area yet minimum energy loss [210]. Another type of structure that is widely under scrutiny is the planar electrode, as utilized in solid-state thin-film batteries [203]. The thin film electrode layer can be viewed as a solid foil, which is a continuous film of active electrode material directly deposited on the current collector. Therefore, the electrolyte can not infiltrate the electrode, rather it is in contact with the planar surface. The lithium species have to travel through this film via diffusion after being deposited at the interface.

The numerical results of the planar electrode and the porous electrode are compared in Figure <ref>. The concentration evolution of the planar electrode shows a clearly defined region of interface with a lithium-rich phase near the separator and lithium-poor phase near the current collector. However, even though more fraction of region near the separator corresponds to the Li-rich phase during the concentration evolution in the porous electrode, no clearly defined interface is detectable along the depth of the electrode. In addition, Figure <ref>(a) shows concentration profiles of the porous electrode with electrode particles aligned perpendicular to the separator, while (b) shows concentration profiles of electrode particles aligned parallel to the separator. The concentration profiles in the latter case oscillate more compared to the former case (red and green curves in Figure <ref>(a) and (b)). This can be justified from the fact that the particles parallel to the separator apparently are stacked layer-wise as shown in Figure <ref>(b1). As the phase separation in a particle starts from the section towards the separator, several particles at the same distance in an affected layer initiate the phase separation simultaneously. Due to the enhanced diffusivity of the electrolyte, the lithium flux traverses through the surrounding of the particles. This initiates the phase separation in the next layer prior to the end of phase separation in the previous layer as shown in Figure <ref>(b3). Therefore, the wave-like pattern can be observed in Figure <ref>(b).

§.§ Comparison with bulk-transport and surface-reaction limited models

[Evolution of SOC and fraction transformed from PFM compared with analytical relations for bulk-transport and surface-reaction limited theories](a) Change in SOC of the electrode and (b) fraction of total particles transformed to Li-rich with time for various particle sizes obtained from the phase-field model (PFM)PFMphase-field model. Here $\tau$$\tau$tortuosity denotes tortuosity, $\rho$$\rho$porosity indicates porosity of the system, and $R$$R$particle whose area is equivalent to the area of the circular particle of radius $R$ denotes particle size whose area is equivalent to the area of the circular particle of radius $R$. Along side, the evolution of SOC for the planar electrode from PFM and the analytical relation for bulk transport limited theory $j_{\textrm{Cot}} \propto t^{-1/2}$ (Cottrell [211] line) are displayed for the comparison in (a). The solid dark-colored lines in (b) are guide to eye for the results from PFM, while light-colored lines are the relations obtained from Johnson-Mehl-Avrami (JMA)JMAJohnson-Mehl-Avrami equation (<ref>).

Simplified kinetics of the system is studied by measuring the net species insertion response of the material under a prescribed concentration at the simulation boundary. The SOC of the porous electrode increases rapidly during the initial phase and shows a plateau subsequently as shown in Figure <ref>(a), which indicates that the flux (rate of change of SOC) of diffusing species of the presented model decays with time. The obtained results are compared with the analytical relations from bulk-transport and surface-reaction limited diffusion [158], which are the standard techniques employed to measure the change in species current in a controlled potential environment in electrochemistry.

[Schematic of bulk diffusion and surface reaction limited transportation dynamics.]Schematic of (a) bulk diffusion and (b) surface reaction limited transportation dynamics. The black arrows indicate dominant driving force, while red arrows illustrate sluggish response. The thick blue line represents electrode surface, while dotted-black line indicates interface between the Li-rich and Li-poor regions.

The bulk-transport limited theory assumes the diffusion of species in the bulk electrode being slower than the deposition at the surface, as illustrated in Figure <ref>(a). Therefore, the transportation rate (or species flux) is controlled by bulk diffusion. A simplistic relation of species flux with time in a bulk-transport limited electrode under a potentiostatic condition is described by Cottrell equation [211],

\begin{equation}
j_{\textrm{Cot}} = k_{\textrm{Cot}} t^{-1/2}, \label{eq:CottrellEquation}
\end{equation}

where $j_{\textrm{Cot}}$$j_{\textrm{Cot}}$Cottrell flux  (in h$^{-1}$ units) is the flux of lithium species measured as the rate of change of SOC and $k_{\textrm{Cot}}$$k_{\textrm{Cot}}$Cottrell flux-time proportionality constant (in h$^{-1/2}$ units) is the proportionality constant associated with operating conditions and electrode properties. The Cottrell equation is commonly considered for planar electrodes in Potentiostatic Intermittent Titration Technique (PITT)PITTPotentiostatic Intermittent Titration Technique to measure the diffusivity [212]. Therefore, the Cottrell equation (<ref>) for $k_{\textrm{Cot}}=\unit[0.75]{h^{-1/2}}$ is plotted with the results obtained from planar electrode for comparison. After the deposition at the particle surface, the species transport through the bulk of the electrode. Due to the lower diffusivity of the electrode, species transport through the bulk of the electrode is the transportation rate-limiting factor, which correlates to the Cottrell equation (<ref>), as shown in Figure <ref>(a). In addition, the results obtained from the porous electrode of various particle sizes are presented. Compared to electrodes of smaller particles, it is evident that the response from the electrode of larger particles ($R=\unit[1.00]{\mu m}$) tends towards the planar electrode and the Cottrell line. Therefore, it can be inferred that the bulk diffusion-limited transportation is more prominent in larger particles compared to the smaller ones.

Another widely recognized theory of the species diffusion is the surface-reaction limited transportation, which considers the bulk diffusion to be more favorable compared to the species deposition at the surface, as shown in Figure <ref>(b). Therefore, the transportation rate is controlled by the surface reaction, in which the nucleation events are a primary focus. Due to dominant bulk diffusion, ideally, the nucleated particles can be regarded as a completely transformed to the Li-rich phase, otherwise consisting of Li-poor phase. An analytical relation of the transformation rate with time can be obtained by invoking few assumptions such as, the nucleation is assumed to occur randomly over the entire untransformed portion and the previously nucleated particle does not influence the likelihood of nucleation around that particle. With these conditions, Johnson-Mehl-Avrami (JMA) equation [213] relates the transformed fraction in two-dimensional systems as,

\begin{equation}
\zeta = 1-e^{N_ct^3}, \label{eq:JMAequation}
\end{equation}

where $\zeta$$\zeta$transformed fraction denotes transformed fraction and $N_c$$N_c$nucleation rate coefficient is the constant associated with nucleation rate. The JMA equation defines the transformation of the particles to follow a sigmoidal shape. Figure <ref>(b) shows results obtained from PFM for three different particle sizes compared with the relations from the JMA equation. As the particle size increases, the transformation rate deviates from the sigmoidal shape considerably. Therefore, the electrode with smaller particles prefers surface reaction limited transportation.

Ultimately, these two theories, utilized in Ref. [158] for intercalation compounds, can be viewed as special cases of the presented numerical model. As an inference from the comparison, the transition of SOC evolution between two different regimes, the bulk-transport and the surface-reaction limited theories can be induced by changing the electrode characteristics such as the size of the particles.

§ INFLUENCE OF MORPHOLOGY DESCRIPTORS ON TRANSPORTATION RATE

§.§ Effect of tortuosity

[Effect of tortuosity $\tau$ on transportation dynamics for different particle sizes $R$]Effect of tortuosity $\tau$ for different particle sizes $R$ as a function of the time required to complete the phase separation process $t_{PE}$, where the SOC stagnates.

The complete process of phase separation is displayed in Figure <ref> for two representative systems: (a) electrode particles aligned perpendicular and (b) parallel to the separator. System (a) consumes $t_{PE}=\unit[10.26]{h}$ for a phase separation process to complete, while system (b) requires $t_{PE}= \unit[23.10]{h}$. As the systems correspond to equal particle sizes and porosity, the distribution of particles can be considered immensely influential on the $t_{PE}$. One of the parameters that quantify the geometrical particle distribution is the tortuosity of the system $\tau$. Based on studies of Refs. [206, 207, 208], the technique to measure the tortuosity is described in Appendix <ref>. For the provided representative cases in Figure <ref>, tortuosity $\tau = 1.28$ attained $t_{PE} = \unit[10.26]{h}$, while $\tau = 2.43$ necessitates $t_{PE} = \unit[23.10]{h}$. Furthermore, for different $\tau$ and maintaining the porosity and the mono-disperse particle sizes of the systems consistent, $t_{PE}$ is  linearly related to tortuosity $\tau$, which is plotted as the red curve in Figure <ref>. Therefore, the increase in tortuosity of the system linearly increases the time required to end the phase separation, $t_{PE}$.

§.§ Effect of particle size

The rate of change of SOC and fraction transformed relate to the particle sizes $R$, as displayed in Figure <ref>(a) and (b) respectively, which shows the numerical result of systems containing several particles of a given size with maintaining equal tortuosity and porosity, where $R$ denotes particle size whose area is equivalent to the area of the circular particle of radius $R$. The system of smaller particles (blue curve) tends towards the surface-reaction limited line in (b), while that of the larger particles (green curve) is situated close to the bulk-transport limited line (black curve) in (a). The discrepancy can be explained by considering the distance diffusing species have to travel after being deposited at the surface of the particle. The diffusing species avail more time to migrate from the surface to reach the center in the larger particles. Therefore, bulk transport is the rate-limiting factor in the larger particles, which is analogous to the bulk-transport limited model. This effect diminishes with the reduction in particle size. Therefore, smaller particles observe higher species flux (rate of change of SOC) in the initial stage and saturate quickly, similar to the surface-reaction limited model.

Figure <ref> shows the effect of different mono-disperse particle sizes as a function of tortuosity $\tau$ and $t_{PE}$ with maintaining porosity. For different particle sizes, as shown in Figure <ref>, the system of smaller particles reaches saturation earlier compared to larger particles. Apart from that, linearity in the $\tau$ and $t_{PE}$ relation is observed for all systems of various particle sizes. Furthermore, $\tau$ vs $t_{PE}$ lines are parallel for different particle sizes with the same porosity. However, the porosity certainly affects the time of $t_{PE}$, which is described in the next section.

§.§ Effect of porosity

[Effect of porosity $\rho$ on transportation dynamics for different tortuosity $\tau$]Effect of porosity $\rho$ for different tortuosity $\tau$ as a function of the time required to complete the phase separation process $t_{PE}$, where the SOC stagnates.

During the formation of electrode packing, the selection of an optimum mass fraction of host material has gathered much attention due to its implications on energy density and performance of the battery. The mass fraction of the host material is directly linked to the porosity of the electrode. Therefore, understanding the effect of variation in porosity on the phase separation dynamics has technological implications. Figure <ref> shows the relation of $\tau$ and $t_{PE}$ with the porosity $\rho$ for a consistent particle size $R$. A linear relation of $t_{PE}$ with $\tau$ can be observed for respective values of porosity.

[Insertion of species in a 3D cathode electrode under a constant concentration boundary condition at the separator]Insertion of species in a cathode electrode with $\unit[4.6]{\mu m}$ depth along with current collector under a constant concentration boundary condition at the separator in three-dimensional space. (a) demonstrates simulation geometry, while others depict phase-separated morphologies for various time-steps. The images correspond to SOC (b) 58% ($t= \unit[0.5]{h}$), (c) 70% ($t= \unit[1.0]{h}$), (d) 85% ($t= \unit[2.6]{h}$), (e) 90% ($t= \unit[3.6]{h}$), and (f) 93% ($t= \unit[4.9]{h}$). The red color represents the lithium-rich phase. The blue color cavity inside the lithium-rich phase corresponds to the lithium-poor phase.

The decrease in porosity results in the increased mass fraction of host material and decreased ion carriers in the form of electrolytes. Hence, the system of lower porosity requires more time $t_{PE}$ to obtain the SOC plateau. Also, the slope of the linear relationship ($\tau$ with $t_{PE}$) is not constant, instead increases with decreasing porosity. This implies that the system of higher $\tau$ exhibits appreciable increment in $t_{PE}$ compared to the lower ones for the same porosity difference.

Finally, Figure <ref> demonstrates the capability of the employed model that the method can be readily extended to simulate three-dimensional systems. The 3D electrode is considered for the numerical experiment to study the transport mechanism. The constant concentration at the separator initiates the formation of the Li-rich phase from the outer surface of the particles, otherwise consists of the Li-poor phase throughout, as shown in (b). The Li-rich phase continues to grow further at the surface of the particles in (c), (d), and (e). Eventually, the lithium-poor phase has been eliminated by the lithium-rich phase in (f). Note that the 3D case shows a significantly different pattern compared 2D case during transportation in the particles. Even though particles are aligned parallel, it shows a transportation mechanism more like the perpendicular 2D case than the parallel 2D case. This is due to the fact that the tortuosity here ($\tau=1.23$) is more like the perpendicular 2D one. Furthermore, the depth of the electrode in 3D is much less than the 2D case in order to minimize the computational efforts. According to the preliminary findings, the depth of the electrode certainly influences the SOC evolution, which will be subject of future investigations.

§ CONCLUSION

A numerical result of the porous electrode is presented to understand the effect of various microstructural properties such as particle size, porosity, and tortuosity on the lithium transport mechanism. In this chapter, the phase separation mechanism in LMO particles has been successfully demonstrated in a multiple particle model system. Ellipsoid-like particles are considered as an example, however, the model can be readily applicable to particles of complicated geometries.  According to the diffusional properties of electrode and electrolyte, a study is conducted on transportation rate dependence with various morphological descriptors of electrode microstructures.

The obtained results suggest that the transportation rate of the system is strongly related to the tortuous pathways formed by the particle orientation, which can be quantified by the tortuosity parameter $\tau$. When controlling the lithium concentration at the separator, the lithium transportation rate is observed to be linearly related to the tortuosity. Furthermore, the slope of this linear relation is independent of the particle size, while the slope alters with a change in the porosity of the electrode. Therefore, the tortuosity, the porosity, and the particle size can be suitable descriptors for the characterization of electrode morphologies.

Furthermore, the evolution of the state of charge of the obtained results for the porous structures of mono-disperse particles are compared with the bulk-transport and surface-reaction limited theories. The results suggest that systems consisting of smaller particles are limited by surface reaction, while larger particles tend towards the bulk-transport limited theory derived for planar electrodes. In order to identify the promising hierarchically structured electrodes, the presented simulation results could be utilized to optimize the experimental efforts.

tocchapterIV Results and Discussion:
Electromigration in metallic conductors

Part IV

Results and Discussion:
Electromigration in metallic conductors

CHAPTER: MOTION OF ISOTROPIC INCLUSIONS

Phase-field methods are extensively applied to study isotropic inclusion propagations under electromigration. In the stimulating work by Mahadevan and Bradley [124], the development of a slit-like feature from an edge perturbation is studied by the phase-field model due to electromigration. The work focused on the physical mechanism of slit growth at the location of a preexisting notch and subsequent propagation transverse to the line. As a result, the change in volume of the inclusion depends strongly on the applied current. Contrarily, experimental evidence also corroborates that the motion of volume-preserving inclusions along the conductor line [214, 16], which instigated interest in modeling community in recent years. Bhate et al. [123] presented the morphological evolution of a geometrical shape under the influence of surface energy, electric and stress fields on single isolated inclusions. Thereafter, Barrett et al. [127] studied the migration, splitting, and coalescence of inclusions along the metallic conductors. Furthermore, Baňas et al. [121] presented results of geometrical shape evolution in three dimensions space. There are instances where propagation of a finger-like slit with shape and volume preserved is observed in SnAgCu solder bumps [177], which has not received attention from a modeling perspective until now. Furthermore, a critical comparison between the solutions obtained from phase-field and sharp-interface methods is lacking.

Therefore, the emphasis of the presented chapter is to compare the results from the phase-field model derived in Chapter <ref>  with the sharp-interface analysis to investigate the selection of slit width and velocity of an initially circular inclusion at an applied electric field. More importantly, the second objective of the study is to investigate the geometrical morphologies of inclusions from the phase-field model, which is otherwise unfeasible with sharp-interface theories. The sharp-interface analysis is considered in section <ref> to derive a critical limit of circular inclusions stability with assuming equal conductivities to the matrix. Afterward, the assumption of equal conductivities is relieved in section <ref>. Thereafter, results from the phase-field numerical study are compared with linear stability analysis in section <ref>. Subsequently, the description of the finger-like slits obtained from the numerical study is presented in section <ref>. Also, the characteristics associated with the slit shapes are derived from sharp-interface model

and a critical comparison with numerical results are also presented.

Finally, the chapter is concluded by its practical implications in section <ref> with subsequent discussion of the important results in section <ref>. Some parts of this chapter are published in the Journal of Electronic Materials [119].

§ SHARP-INTERFACE ANALYSIS FOR CIRCULAR ISLANDS

[Schematic of a circular island, subjected to an external electric field $E_\infty$, in an infinite conductor]Schematic of a circular island, subjected to an external electric field $E_\infty$, in an infinite conductor domain.

The stabilty analysis of a circular dislocation loop propagation, under the assumption of a steady-state, is derived by Yang et. al [162]. For brevity, the basic steps of the derivations are casted here for islands migrating in a conductor. In the present context, the term, island describes a cluster of material entrapped in the metallic conductor consists of equivalent conductivity to that of the conductor.

Consider an island, translating along the length of the conductor under an external electric field, as shown in Figure <ref>. The migration of the island is due to the result of a mass transport flux, which is induced by surface electromigration, and by the capillarity given by the Nernst-Einstein relation:

\begin{equation}
J_s = \frac{D_s \delta_s}{\Omega k_B T} \Big( -eZ_sE_t + \Omega \gamma_s \frac{\textrm{d} \kappa_s}{\textrm{d} s} \Big), \label{eq:EM1SharpInterfaceFlux}
\end{equation}

where $J_s$$J_s$surface atomic flux represents the number of species passing per length of island surface per unit time, $D_s$ is the surface diffusivity, $\delta_s$ is the thickness of the surface layer, $\Omega$ denotes the atomic volume, $e$ is the electron charge, $Z_s$$Z_s$effective valence represents the effective valence at the surface, $k_B$ is the Boltzmann constant, $T$ denotes the absolute temperature, $\gamma_s$ is the surface energy, $\kappa_s$ represents the local curvature of the surface, and $s$$s$arc length along the island/inclusion surface denotes the arc length along the island surface.
The capillary-mediated flux is driven by the gradients of curvature, along the island surface, while the electromigration flux is dictated by the tangential component of the electric field along the island surface $E_t$$E_t$tangential component of the electric field along the island/inclusion surface, given by

\begin{equation}
E_t = - E_{\infty} \textrm{ sin } \theta, \label{eq:EM1sharpInterfaceElectricalPotentialIsland}
\end{equation}

which is responsible for the island transmission [215]. $E_{\infty}$$E_{\infty}$applied electric field is the applied electric field, $\theta$ is the tangent angle of the island surface.

Mass conservation relates the surface divergence of the flux to the normal velocity$V_n$velocity along the surface normal as follows:

\begin{equation}
\frac{\textrm{d}J_s}{\textrm{d}s} = \frac{V_n}{\Omega}. \label{eq:EM1MassConservationNormalVelocity}
\end{equation}

By assuming a steady state, the normal velocity can be related to the velocity$V$velocity along the external electric field along applied electric field as

\begin{equation}
V_n = V \textrm{ cos } \theta. \label{eq:EM1normalVelocity}
\end{equation}

The $\kappa_s$curvature along the island/inclusion surfacecurvature along the island surface is the fraction of change in tangent angle to the length of line segment, which can be expressed as,

\begin{equation}
\kappa_s = \frac{\textrm{d}\theta}{\textrm{d}s}.\label{eq:EM1surfaceCurvatureIsland}
\end{equation}

Substituting Eqs. (<ref>), (<ref>), (<ref>), and (<ref>) into the mass conservation equation (<ref>), the resultant expression can be written as,

\begin{equation}
V\textrm{ cos } \theta  = \frac{D_s \delta_s}{ k_B T} \left( e Z_s E_\infty \textrm{ cos } \theta \kappa_s + \Omega \gamma_s \kappa_s \frac{\textrm{d} }{\textrm{d} \theta} \left( \kappa_s \frac{\textrm{d}  \kappa_s}{\textrm{d} \theta} \right) \right). \label{eq:EM1resultantExpression}
\end{equation}

This equation can be non-dimensionalized using the relation,

\begin{equation}
\left[1/\kappa_s \right] = \left[ 1/\hat{\kappa}_s \right] \frac{D_s \delta_s}{k_B T V} eZ_s E_\infty \label{eq:EM1nonDimensional}
\end{equation}

Here $\hat{\kappa}_s$$\hat{\kappa}_s$non-dimensional curvature of the island surface denote the non-dimensional curvature of the island. Substituting Eq. (<ref>) into Eq. (<ref>), the shape equation can be expressed in a non-dimensional form as,

\begin{equation}
\frac{\textrm{d} }{\textrm{d} \theta} \left( \hat{\kappa}_s \frac{\textrm{d}  \hat{\kappa}_s}{\textrm{d} \theta} \right) + \chi \left( 1 - \frac{1}{\hat{\kappa_s}} \right) \textrm{ cos } \theta =0. \label{eq:EM1differentShapes}
\end{equation}

where the dimensionless number $\chi$$\chi$dimensionless number associated with the steady state shapes of an island is defined as,

\begin{equation}
\chi = \left( \frac{D_s \delta_s}{k_B T V} eZ_s E_\infty \right)^2 \frac{eZ_sE_\infty}{\Omega \gamma_s}. \label{eq:EM1ChiDifferentShapes}
\end{equation}

The non-linear ordinary differential equation (<ref>) determines the steady-state shapes. Note that, heretofore, the expressions are equally applicable to non-circular islands. Moving on to circular shaped islands, their stability can be analyzed by introducing perturbations to the surface in the form,

\begin{equation}
\hat{\kappa}_s = 1 + \varepsilon \varkappa (\theta) \label{eq:EM1islandPerturbation}
\end{equation}

where $\varepsilon$$\varepsilon$perturbation parameter is the small perturbation parameter and $\varkappa (\theta)$$\varkappa$perturbation function denotes the perturbation function. Substituting Eq. (<ref>) into Eq. (<ref>), the growth of perturbations can be expressed as,

\begin{equation}
\frac{\textrm{d}^2 \varkappa}{\textrm{d}\theta^2} + \chi \varkappa \textrm{ cos } \theta =0. \label{eq:EM1MathieuEquation}
\end{equation}

This non-trivial expression is similar to Mathieu equation. A MATLAB based subroutine is developed with bpv4c  algorithm to seek eigenvalue of Eq. (<ref>). The second order ordinary differential equation splits into two first order equations with initial guess function and boundary conditions are employed in the form of a cosine function. The lowest eigenvalue is estimated as $\chi=10.65$.

If the island does not collapse, and keeps its circular shape of radius $R_i$$R_i$radius of circular island/inclusion, the steady-state velocity$V_0$steady-state velocity of an island/inclusion [74] can be determined by using $\kappa_s =1/R_i$ and $\textrm{d} \kappa_s / \textrm{d}s =0$. Therefore, substituting these relations in Eq. (<ref>), the velocity of a circular island can be expressed as,

\begin{equation}
V_0 = - \frac{D_s \delta_s}{R_ik_BT} eZ_s  E_{\infty} \label{eq:EM1steadyStateCircularVelocity}
\end{equation}

Therefore, the dimensionless number in Eq. (<ref>) can be expressed for steady-state circular islands of the form,

\begin{equation}
\chi_0 = \frac{eZ_sE_\infty R_i^2}{\Omega \gamma_s}. \label{eq:EM1chi0}
\end{equation}

This equation represent relative strength of the capillarity $\Omega \gamma_s/R_i^2$ and the electromigration force $eZ_s E_\infty$. The capillarity attempts to reduce any curvature gradient in an ultimate aim to form a circular shape, while the atomic transfer due to electromigration in the direction of electron wind alters the preexisting circular shape of the island. Therefore, the island shape is controlled by the parameter $\chi_0$$\chi_0$dimensionless number associated with the stability of circular island/inclusion, which represents the competition between the electromigration force and the capillarity. A critical limit of the electric field, at which a circular island maintains its shape is defined by $\chi_{0\textrm{ct}}=$$\chi_{0\textrm{ct}}$critical value of $\chi_0$, above which circular island/inclusion is unstable 10.65. Above the critical value, i.e., when the electron wind force prevails over the capillary force, a circular island is unstable and undergoes shape bifurcation.

Heretofore, a study of the homogeneous inclusion is considered. Meaning, the conductivity of the circular inclusion is equivalent to the metallic conductor. To differentiate, the homogeneous inclusions are termed as the islands. However, the presented theory needs to be extended in order to capture the behavior of insulating voids and different conductivity material entrapped into metallic conductors, which has technological relevance. Therefore, in the next section, a theory is developed to determine the stability limit of circular heterogeneous inclusions.

§ LINEAR STABILITY ANALYSIS OF CIRCULAR INCLUSION

A direct projection of external electric field at the infinity is considered on the surface of inclusion in previous section <ref>(see Eq. (<ref>)). This can be justified for homogeneous inclusions. However, for heterogeneous inclusions, a general expression of electrical potential in a matrix containing a circular inclusion is determined. Then, this expression is extended to a circular inclusion with small perturbations on the surface. Finally, a stability limit of a circular inclusion is determined from the basic species transport relation. The presented theory is in coherence with the work of Ref. [112]. For the brevity and completeness, the approach with necessary assumptions is described.

§.§ General expression of electrical potential in a matrix containing a circular inclusion

Consider a circular inclusion of radius $R_i$ in an infinite matrix subject to an external electric field, $E_{\infty}$. In the present study, the plane polar coordinate system $(r, \vartheta)$$r, \vartheta$polar coordinates of a point is considered with the origin at the center of the inclusion. $r$ represents a distance of a point in the system (i.e., matrix + inclusion) and $\vartheta$ corresponds to an angle formed by the point under consideration at the origin with the direction of electric field. The electrical potential at any far point from the inclusion is stated as,

\begin{equation} \label{eq:linearStabilityElectricPotential}
\phi (r,\vartheta) \rvert_{r \rightarrow \infty} = - E_{\infty} r \textrm{ cos } \vartheta.
\end{equation}

It is easy to perceive that the inclusion of conductivity $\sigma_{\textrm{icl}}$ disrupts the electrical potential distribution in the matrix of different conductivity, $\sigma_{\textrm{mat}}$, inside and nearby the inclusion. The spatial distribution of the electrical potential can be obtained by Laplace's equation, in cylindrical coordinates with azimuthal symmetry,

\begin{equation}
\frac{1}{r}\frac{\partial}{\partial r} \left( r \frac{\partial \phi (r,\vartheta)}{\partial r}\right) + \frac{1}{r^2}  \frac{\partial }{\partial \vartheta}\left( \frac{\partial \phi (r,\vartheta)}{\partial \vartheta} \right) = 0.
\end{equation}

The general solution of the above equation can be derived by the variable separation method,

\begin{equation}
\phi (r,\vartheta) = \sum^{\infty}_{l=1} \left( A_l r^l +\frac{B_l}{r^l} \right) (C_l \textrm{ cos }l\vartheta + D_l \textrm{ sin }l\vartheta).
\end{equation}

$A_l, B_l, C_l,$ and $D_l$ are the real coefficients, which need to be determined by known conditions. At the far distance, above equation of the electrical potential should correspond to Eq. (<ref>). Thus the equation reduced to the form,

\begin{equation}
\phi(r,\vartheta) = A r\textrm{ cos } \vartheta + \frac{B }{r} \textrm{ cos } \vartheta,
\end{equation}

containing only cos $\vartheta$ term. Furthermore, as the origin is at the center of the inclusion, the second term diverges at the origin. Therefore, electrical potential outside and inside of the sphere is expressed separately as,

\begin{eqnarray}
\phi_{\textrm{in}}(r,\vartheta) &=& A_{\textrm{in}} r \textrm{ cos } \vartheta,\\
\phi_{\textrm{out}}(r,\vartheta) &=& A_{\textrm{out}} r \textrm{ cos } \vartheta +  \frac{B_{\textrm{out}}} {r} \textrm{ cos } \vartheta.
\end{eqnarray}

Comparing above equation with the asymptotic limit at the far distances, the expression $A_{\textrm{out}} = - E_{\infty}$ can be obtained. It is easy to notice that, from the continuity equation, on the inclusion surface,

\begin{equation} \label{eq:linearStabilitycontinuityConditions}
{ \left. \frac{\partial \phi_{\textrm{in}}}{\partial r}\right|}_{r=R_i} = \beta
{ \left. \frac{\partial \phi_{\textrm{out}}}{\partial r}\right|}_{r=R_i} \hspace{0.5cm} \textrm{and} \hspace{0.5cm} { \left. \frac{\partial \phi_{\textrm{in}}}{\partial \vartheta}\right|}_{r=R_i} =
{ \left. \frac{\partial \phi_{\textrm{out}}}{\partial \vartheta}\right|}_{r=R_i} .
\end{equation}

Here $\beta = \sigma_{\textrm{mat}}/\sigma_{\textrm{icl}}$$\beta$conductivity ratio, $\sigma_{\textrm{mat}}/\sigma_{\textrm{icl}}$. The resultant electrical potential reads,

\begin{equation}
\phi(r,\vartheta) =     \left\{ \begin{array}{ccl}
\frac{-2\beta}{(1+\beta) }E_{\infty} r \textrm{ cos }\vartheta, & \textrm{  for}& r \leq R_i, \\
-E_{\infty} r \textrm{ cos } \vartheta + \frac{E_{\infty} R_i^2}{r} \frac{(1-\beta)}{(1+\beta)} \textrm{ cos } \vartheta & \textrm{  for} & r \geq R_i.
\end{array}\right. \label{eq:linearStabilityflux}
\end{equation}

Note that both expressions are satisfied at the surface of the inclusion.

§.§ Electrical potential in a matrix containing a small perturbed circular inclusion

Using the conformal mapping to express complex potential (Eq. (<ref>)) for the circular inclusion at point $z = x + iy$$i$complex number in complex form,

\begin{eqnarray}
\Phi_{in}(z)&=& \frac{-2\beta }{(1+\beta)} E_{\infty} z,\\
\Phi_{out}(z)&=& - E_{\infty} z + \frac{(1-\beta)}{(1+\beta)} E_{\infty} \frac{R_i^2}{z}.
\end{eqnarray}

Here $\Phi(z) = \phi + i \psi$, where $\psi$$\psi$orthogonal function to the electrical potential is the orthogonal function to the electrical potential. Heretofore, the linear stability of a perfectly rounded inclusion is considered. Now, a small perturbation $(\varepsilon \ll 1)$ is introduced to the surface of inclusion. In the present part, an expression of electric field is derived for the matrix containing inclusion of circular shape with small perturbations around the surface. The mapping of inclusion surface,

\begin{equation}
z(R_i, \vartheta) = x + iy = R_i(\textrm{cos }\vartheta + i \textrm{ sin }\vartheta ) + \varepsilon \sum^{\infty}_{n=0} a_n (\textrm{cos } n \vartheta - i \textrm{ sin }n \vartheta ),
\end{equation}

here $a_n$ are real constants$a_n$surface perturbation coefficient corresponds to n$^{\textrm{th}}$-term, which associated with the strength of surface perturbations. The expression for the electrical potential needs to be updated, as the small perturbations on the inclusion surface modify the electrical potential distribution. From the arguments based on a convergence of the electrical potential inside and outside of the inclusion, the expression of electrical potential can be given by,

\begin{eqnarray}
\Phi_{in}(z)&=& \frac{-2\beta }{(1+\beta)} E_{\infty} z + \varepsilon \sum^{\infty}_{n=0} C_n z^n,\\
\Phi_{out}(z)&=& - E_{\infty} z + \frac{(1-\beta)}{(1+\beta)} E_{\infty} \frac{R_i^2}{z} + \varepsilon \sum^{\infty}_{n=0} D_n z^{-n}.
\end{eqnarray}

Here $C_n$ and $D_n$ are coefficients, which are determined by substituting $z$ in the above equation and applying continuity conditions, Eq. (<ref>), on the surface of the inclusion, after transforming, from $\partial \phi_{in}/\partial r = \beta (\partial \phi_{out}/\partial r) $ to $\psi_{in} = \beta \psi_{out}$ by the Cauchy-Riemann condition (i.e., $\partial \phi / \partial r= \partial \psi / \partial \vartheta$) and from $\partial \phi_{in}/\partial \vartheta = \partial \phi_{out}/\partial \vartheta $ to $\phi_{in} = \phi_{out}$. Therefore, the resultant electrical potential on the surface of inclusion,

\begin{eqnarray}
\Phi_{in} &=& \Bigg{\lbrace}-\frac{2\beta}{(1+\beta)} E_{\infty} R_i \, e^{i\vartheta}  - \varepsilon \frac{2\beta}{(1+\beta)} E_{\infty} \sum^{\infty}_{n=0} a_n e^{-in\vartheta}\\
&&\hspace{5.5cm} - \varepsilon \frac{2 \beta(1-\beta)}{(1+\beta)^2}  E_{\infty} \sum^{\infty}_{n=0} a_n e^{in\vartheta}\Bigg{\rbrace}, \nonumber\\
\Phi_{out} &=& -E_{\infty} R_i \, e^{i \vartheta} + \frac{(1-\beta)}{(1+\beta)} E_{\infty} R_i \, e^{-i \vartheta} - \varepsilon \frac{4 \beta}{(1+\beta)^2}E_{\infty}\sum^{\infty}_{n=0} a_n e^{-in\vartheta}.\label{eq:linearStabilityelectricalPotential}
\end{eqnarray}

Here $e^{i\vartheta} = \textrm{ cos }\vartheta +i \textrm{ sin }\vartheta$.

§.§ Linear stability analysis of transient inclusion

The mapping function $z$ can be expressed for a transient inclusion, which drifts in the matrix. The mapping function of the inclusion, whose position changes with time can be given by,

\begin{equation}\label{eq:linearStabilitysurfaceMapping}
z(R_i,\vartheta, t) = R_i\, e^{i\vartheta} + S(t) + \varepsilon \sum^{\infty}_{n=0} a_n (t) \, e^{-in\vartheta}.
\end{equation}

where $S(t)$ is the spatial displacement function related to path traveled by the inclusion$S$displacement function. In the present case, the coefficients $a_n$ vary with time as the inclusion shape transforms as it evolves.

[Schematic of a line segment of a surface of inclusion.]Schematic of a line segment of a surface of inclusion.

Consider a small segment on the inclusion surface as shown in Figure <ref>. Length of the segment can be expressed as,

\begin{equation}\label{eq:linearStabilitysegmentLength}
\textrm{d}s = \vert \textrm{d}z\vert =\left\vert \frac{\textrm{d}z}{\textrm{d} \vartheta} \frac{\textrm{d}\bar{z}}{d\vartheta} \right\vert^{1/2} d\vartheta = R_i\left( 1 -   \frac{\varepsilon}{R_i} \sum^{\infty}_{n=0} na_n \textrm{ cos }(n+1) \vartheta \right) d\vartheta.
\end{equation}

Here, $\bar{\bullet}$ is a complex conjugate of the respective function. The curvature of the segment can be expressed as,

\begin{equation}
\kappa_s = \frac{\textrm{d} \theta}{ \textrm{d} s} = \frac{\textrm{d}}{\textrm{d}s}\left(\textrm {tan}^{-1} \left(- \frac{\textrm{d}x/\textrm{d}s}{\textrm{d}y/\textrm{d}s} \right) \right).
\end{equation}

where $\theta$ is an angle formed by line segment with the perpendicular to external electric field. Simplifying above equation by substituting the expressions for $\textrm{d}y/\textrm{d}s$, $\textrm{d}x/\textrm{d}s$ and $\textrm{d}s/\textrm{d}\vartheta$ from the Eqs. (<ref>) and (<ref>), one obtains

\begin{equation}\label{eq:linearStabilitycurvatureMapping}
\kappa_s = \frac{1}{R_i} \left( 1 + \frac{\varepsilon}{R_i} \sum^{\infty}_{n=0} n(n+2) a_n \textrm{ cos }(n+1)\vartheta \right).
\end{equation}

Now, the electric field can be evaluated from Eqs. (<ref>) and (<ref>) as,

\begin{equation}\label{eq:linearStabilityElectricFieldOnPerturbedVoid}
\begin{aligned}
E_t = - \frac{\textrm{d} \phi_{out}}{\textrm{d}s} = \Bigg\lbrace- \frac{2\beta}{(1+\beta)}&E_{\infty} \textrm{ sin } \vartheta - \varepsilon \frac{\beta (3-\beta)}{(1+\beta)^2} \frac{E_{\infty}}{R_i} \sum^{\infty}_{n=0} n a_n \textrm{ sin }n \vartheta \\
&- \varepsilon \frac{\beta}{(1+\beta)} \frac{E_{\infty}}{R_i} \sum^{\infty}_{n=0} a_{n-2} (n-2) \textrm{ sin } n\vartheta\Bigg\rbrace,
\end{aligned}
\end{equation}

It is easy to verify that by putting $\varepsilon=0$ in the above equation, the expression of the electric field on a perfectly circular inclusion is obtained, which is widely accepted form in the sharp-interface analysis [119, 110].

The normal velocity at any point on the line segment is given by,

\begin{equation}
V_n = \Omega \frac{\textrm{d}J_n}{\textrm{d}s} = V_2 \textrm{ cos }(270^\circ + \theta) +V_1 \textrm{ sin }(270^\circ + \theta) = \frac{\textrm{d}y}{\textrm{d}t} \frac{\textrm{d}x}{\textrm{d}s} - \frac{\textrm{d}x}{\textrm{d}t} \frac{\textrm{d}y}{\textrm{d}s}
\end{equation}

Substituting Eqs. (<ref>) and (<ref>) in the above equation and integrating, the resultant equation can be given by,

\begin{equation}\label{eq:linearStabilityOmegaJ}
\Omega J_s = - \dot{S}R_i\textrm{ sin }\vartheta - \varepsilon R_i \sum^{\infty}_{n=0} \frac{\dot{a_n}}{(n+1)} \textrm{ sin }(n+1)\vartheta + \varepsilon \dot{S} \sum^{\infty}_{n=0} a_n \textrm{ sin } n\vartheta .
\end{equation}

Substituting Eqs. (<ref>), (<ref>), and (<ref>) into the basic mass transport equation,

\begin{equation}
\Omega J_s = \frac{D_s \delta_s}{k_B T}\left( -eZ_s E_t +\Omega \gamma_s \frac{\textrm{d}\kappa_s}{\textrm{d}s} \right).
\end{equation}

In addition, the obtained equation is equivalently applicable to the perfectly circular inclusion. Therefore, by substituting $\varepsilon=0$, one finds a relation of a velocity of a perfectly circular inclusion,

\begin{equation}
\dot{S} = V_0= -\frac{D_s \delta_s}{k_B T} \frac{2\beta}{(1+\beta)}\frac{eZ_sE_{\infty}}{R_i}. \label{eq:linearStabilityinclusionVelocity}
\end{equation}

and the equation of evolution of inclusion surface reduced to,

\begin{equation}
\begin{aligned}
-\dot{a}_{n-1}= \frac{D_s \delta_s}{k_BT} \frac{\beta}{(1+\beta)} \frac{eZ_sE_{\infty}}{R_i^2}n \Bigg( (n-2) a_{n-2} &-\frac{(1+\beta)}{\chi_0 \beta} n(n^2-1) a_{n-1}\\
&+ \left(\frac{(3-\beta)}{(1+\beta)}n+2 \right)a_n \Bigg),
\end{aligned} \label{eq:linearStabilityChi}
\end{equation}

where, a dimensionless parameter is defined as $\chi_0 = eZ_sE_{\infty}R_i^2/(\Omega \gamma_s).$ At the equilibrium of inclusion shape, i.e. $\dot{a}_{n-1}=0$, a set of an infinite number of homogeneous linear equations are obtained of the form,

[Convergence of the numerical solution for linear stability analysis]Convergence of the numerical solution for linear stability analysis. $N$ indicates the number of equations considered for the study. The error is calculated on the basis of relation obtained from 36 equations.

\begin{equation}
(n-2)a_{n-2} - \frac{(1+\beta)}{\beta \chi_{0ct}} n (n^2-1)a_{n-1} + \left( n\frac{3-\beta}{1+\beta} +2 \right) a_n =0,\hspace{1cm} n\geq2.\label{eq:linearStabilitySharpInterface}
\end{equation}

For all integer values of $n\geq2$, a homogeneous linear equation can be obtained in the form of coefficients, $a_{n-2}$, $a_{n-1}$, $a_{n}$, $\beta$, and $\chi_{0ct}$. The purpose of these equations is to find a relation between $\beta$ and $\chi_{0ct}$ for non-trivial coefficients, $a_n$'s. Only a subset of an infinite number of equations is considered due to restrictions on the computational and analytical approaches. In other words, finite set of equations are considered with assumption that coefficients beyond $a_N$ (i.e., $a_{N}$, $a_{N+1}$,...) are zero, where $N$ indicates the number of equations considered for the study. The analytical solution can be obtained by substitution of variables up to five equations. Beyond that, the difficulty to obtain analytical relation increases exponentially with an increase in the number of equations. Therefore, a numerical tool is developed on MATLAB platform by using vpasolve algorithm. The convergence study of the resultant relation is plotted in the Figure <ref>. The derived relation between $\beta$ and $\chi_{0ct}$ in the present section (Eq. (<ref>)) is considered further for the comparison with the results obtained from the phase-field method in Figure <ref>.

§ PHASE-FIELD RESULTS OF INITIALLY-CIRCULAR INCLUSION MIGRATING UNDER ISOTROPIC DIFFUSION

[Dependence of an inclusion equilibrium shape on the dimensionless parameter, $\chi_0$ and the conductivity ratio $\beta$]Dependence of an inclusion equilibrium shape on the dimensionless parameter, $\chi_0$ and the conductivity ratio $\beta$. The black solid curve corresponds to $\chi_{0ct}$, which is obtained by linear stability analysis. This line identifies whether an inclusion to retain its circular shape or evolves to subsequent collapse. The black dashed line indicates lines of action on which the phase-field simulations are obtained. The inset of images are the snapshots from the phase-field model.

The competition between the EM wind force and the capillary force, defined by the non-dimensional parameter $\chi_0$, determines the stability of the circular inclusion. Accordingly, there exists a critical value $\chi_{0ct}$, above which the circular inclusion can no longer be stable and loose its circular shape. The theoretical value of the dimensionless parameter, $\chi_{0ct}=10.65$ , for homogeneous island ($\beta=1$) is obtained by computing the eigenvalue of Mathieu-type equation (<ref>), obtained in section <ref>. The derived value is in agreement with linear stability analysis presented in section <ref>. Furthermore, as shown in Figure <ref>, $\chi_{0ct}$ shows a strong dependency with conductivity ratio, $\beta$.

Even though, this theory provides information regarding the stability of circular inclusion. However, it serves no information about the shape of the inclusion after the collapse from the circular shape. These inferences can be understood by the analysis of the phase-field simulations.

Inset of images from phase-field simulations in Figure <ref> shows the influence of the non-dimensional parameters $\chi_0$ and $\beta$ on inclusion equilibrium shapes. On the one hand, the competition between electromigration wind force and the capillary force decides the consequences on the inclusion shape. The non-dimensional parameter $\chi_0$ determines the relative strength of the electromigration force compared to the capillarity. On the other hand, $\beta$ identifies maximum and minimum values of $E_x$$E_x, E_y$local electric field components in x and y-directions and $E_y$ for a known external electric field $E_{\infty}$, which is responsible for the change of local electric field distribution. The significance of these two parameters is discussed in the following paragraphs. The simulation parameter set considered for this study is summarized in Table <ref>.

Values of parameters considered for the simulation sets of an isotropic inclusions to study significance of $\chi_0$ and $\beta$.

$\chi_0$   $\beta$

Set: 1   $5.0$ to $60.0$   $3$

Set: 2   $7.5$   $1$ to $10000$

Set: 3   $15.0$   $1$ to $10000$

§.§ Significance of $\chi_0$

[Representative case of steady-state circular inclusion and circular inclusion to slit transition](a) A circular inclusion migrates by preserving its shape, and (b) evolves its shape to the finger-like slit. The results correspond to the dimensionless parameter $\chi_0=5.39$ for (a), and $\chi_0=26.99$ for (b) keeping the conductivity contrast constant, $\beta=3$. The arrows show the direction of the electron wind, from the cathode to the anode. Here $u$ is the half-slit width$u$half-slit width and $w$ denotes the line width of the conductor$w$line width of the conductor.

The shape of the inclusion is governed according to the competition between the electromigration wind force and the capillary force, as it is evident from the non-dimensional parameter $\chi_0$.
The capillary force seeks a uniform curvature and consequently diffuses species to reduce any curvature gradient along the inclusion surface.
Electromigration, on the other hand, instigates an atomic transport in the direction of the electron wind force, leading to its shape alterations.
In the present case, this corresponds to the diffusion of the species, from the right towards the left surface.

The migration of the circular inclusion, at a low electromigration force, is shown in Figure <ref>(a).
From the phase-field simulations, the inclusion surface is extracted as an isoline of value $c=0.5$.
The diffusion of the atomic species, in the direction of the wind force, induces the drift of the inclusion, in the direction of the electric field. Inherently, the capillary force suppresses the presence of any curvature gradient along the inclusion surface. The dominant capillary force, however, encourages the inclusion to preserve the circular shape during the course of the drift as shown in Figure <ref>(a). Any non-uniformity of the curvature is almost instantaneously eliminated by the capillarity. Therefore, the circular inclusion is the most preferred shape by the capillarity.

A dominant electromigration force, on the other hand, leads to a biased mass transport towards the anode side, thereby producing a sharp protrusion on the cathode side of the inclusion (see Figure <ref>(b)).
The differential curvature induces a counter mass transport towards the cathode end of the inclusion.
The curvature-mediated mass transport is evident in the shortening rear end of the inclusion.
However, owing to the significantly higher electromigration force, the protrusion on the inclusion surface, at the cathode end, increases in length, and ultimately forms a narrow finger-like slit.

§.§ Significance of $\beta$ and the current crowding

[Velocity of the centroid of the inclusions along the direction of external electric field, $V_0$ for conductivity ratios, $\beta=1, 3,$ and 10000]Velocity of the centroid of the inclusions along the direction of external electric field, $V_0$ for conductivity ratios, $\beta=1, 3,$ and 10000, and the dimensionless parameter $\chi_0=7.50$.  The inset of images are the snapshots from the phase-field model in equilibrium.

To understand the significance of $\beta$, simulations along two horizontal dotted lines shown in Figure <ref> are considered. According to linear stability analysis, all inclusions on the lower line should migrate maintaining their circular shapes, while shape alterations are expected on the right part of the upper line.

From the phase-field simulations, the migration of the initially-circular inclusions for lower electromigration force ($\chi_0=7.5$) are shown in Figure <ref>. Due to capillarity mediated diffusion, circular inclusions are unchanged from its initial configuration irrespective of conductivity ratio $\beta$. However, the velocities of the inclusions are modified by $\beta$, which is in agreement with the literature [119, 110] and linear stability analysis (see Eq. (<ref>)). In fact, the adjustment of the inclusion velocity approximately corresponds to the factor of $2\beta/(1+\beta)$. This infers that the $\beta$ increases velocity of the circular inclusions, as inclusions of heterogeneous conductivities ($\beta>1$) travels faster compared to its homogeneous ($\beta=1$) counterpart. Therefore, the increase in conductivity ratio improves propagation speed, up to two times for a circular inclusion.

[Change in local electric field distribution due to current crowding near inclusion surface]The change in local electric field distribution during the evolution of the inclusion. The distribution of electric field components $E_x$ in first column and $E_y$ in second column are shown for $\beta=1, 3,$ and $10000$ from top to bottom respectively. The electric field components are normalized by the external electric field $E_{\infty}$. The electric field distribution remains unchanged for $\beta=1$. On the contrarily, although it remains nearly constant at the distance from the inclusion, while inside and around the inclusion, electric field changes effectively for $\beta=3$ and $10000$. The contour of maximum value of $E_x$(=1.5) is divided into two regions inside the inclusion for $\beta=3$. In addition, $E_x$ enhances local electric field up to two times, while $E_y$ forms a quadrupole pattern around the inclusion of $\beta=10000$.

During the inclusion migration at higher electromigration force ($\chi_0=15$), when electron wind passes from the periphery, it interacts with the surface. Therefore, the electromigration force perceives the length of the inclusion along the traverse direction to the electric field as a resistance in the path of electron wind. The higher electromigration force reduces the resistance in the path by altering inclusion shape to the formation of a finger-like slit (see inset at the right side of Figure. <ref>). However, these arguments are only valid for the homogeneous systems, meaning, of similar conductivities ($\beta=1$), because of uniform electric field distribution as shown in Figure <ref>(a1) and (a2). For heterogeneous systems due to different conductivities in the matrix and the inclusion, the distribution of electric field near the inclusion disrupted, which significantly contributes to its shape. Figure <ref>(b) and (c) show the local change in electric field distribution in the vicinity of the inclusion with help of the electric field components, $E_x$ and $E_y$, in the longitudinal and traverse directions respectively.