Effect of tortuosity, porosity, and particle size on phase-separation dynamics of ellipsoid-like particles of porous electrodes: Cahn–Hilliard-type phase-field simulations

Improvements concerning the capacity and rate-capability of battery systems can not only be achieved by choosing suitable materials, but also by tailoring the electrode morphologies. Thus, a simulation study is performed to understand the influence of various microstructural properties such as particle size, porosity, and tortuosity on the transport mechanism. In this work, the classical Cahn–Hilliard model is extended to a multiple particle model system. We consider ellipsoid-like particles 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 the electrode structures. Under Dirichlet conditions for concentration, simulation results predict a linear dependence of the characteristic time on tortuosity. These lines are converging with variation in particle size at higher tortuosity values, while they are diverging with variation in porosity. Furthermore, the results suggest that systems consisting of smaller particles are limited by surface reaction while larger particles tend toward the bulk-transport limited theory derived for planar electrodes.

Besides, Orvananos et al [40] demonstrated the impact of electrode architecture in terms of ionic and electronic connectivities by considering interactions between active particles. Furthermore, Vasileiadis et al [41] investigated that the tortuosity, particle sizes, porosity, and electrode thickness influence the capacity of Li 4 Ti 5 O 12 (LTO) electrodes as a function of C-rate. However, most of the studies utilize an approximation by the Bruggeman relation [42] for the calculation of tortuosity, instead of estimating the tortuosity based on the electrode microstructure under investigation [43][44][45]. In addition, these porous electrode PFMs are focused on LiFePO 4 or Li 4 Ti 5 O 12 electrode materials. Relatively few models have been developed for LMO materials [22,25,29,30], which are limited to single particles. Therefore, this article 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 [43].
In the present work, we extend the Cahn-Hilliard-type PFM to characterize the transportation rates of porous electrode structures as shown in figure 1 with application of two-phase coexistence in LMO particles. In this initial effort, neglecting the reaction kinetics at the electrode-electrolyte interface, we focus on the diffusional properties of the electrode and the electrolyte. In addition, we avoid an explicit treatment of mechanical effects due to misfit strains at the phase boundary. Instead, the emphasis is placed on the influence of various microstructural properties such as particle sizes, porosity, and tortuosity of electrode consisting of ellipsoid-like particles with a focus on the transport mechanism. The present article is organized as follows: section 2 illustrates the employed numerical method based on Cahn-Hilliard equation for a porous electrode. Few typical cases of insertion dynamics are discussed with highlighting the underlying mechanism in section 3. Also, the obtained results are analyzed and compared with known relations. Furthermore, the article is concluded by a brief discussion on the applicability of the presented results in section 4.

Model formulation
A concentration parameter c (x, t) ∈ [0, c max ] is introduced to indicate the local state of the concentration with time t and spatial coordinates x, where c max denotes maximal concentration of lithium species in the electrode material. The process of non-equilibrium thermodynamics, the lithium insertion and the extraction, is governed by the system free energy of the form where f chem CH is the chemical free energy density, f grad CH denotes the gradient energy density, Ω represents the volume of the simulation domain.
The model focuses on a representative elementary volume of an electrode. The set of indicator parameters ψ = [ψ 1 , . . . , ψ N ], ψ a ∈ {0, 1}, ∀a ∈ {1, . . . , N} is defined to recognize different phases. For instance, ψ 1 = 1 (and ψ a = 0, where a = 1) represents electrode particles, while ψ 2 = 1 (and ψ a = 0, where a = 2) indicates the electrolyte. Thus, the indicator parameter identifies multiple regions of the same physical properties as a single phase. The chemical-free energy density f chem CH is a function of bulk properties, such as local concentration, temperature, and orientations. In general, the chemical free energy density can be expressed as where h(ψ a ) is the interpolation function and f chem a denotes the chemical-free energy density of the phase a. Here, as the phases are defined by a step function of indicator parameters, meaning, the interface between the different phases is sharply defined, a first-order interpolation h(ψ a ) = ψ a can be utilized [46]. At every position in the simulation domain Σψ a = 1 and Σh(ψ a ) = 1, while h(ψ a ) can either be 0 or 1 depends on the value of ψ a .
The coexistence of lithium-rich and poor phases in several electrode compounds is of scientific importance to study the phase separating behavior. To capture this effect, a regular solution model is considered in the form of where α a and α a are the regular solution constants related to material properties, k B is the Boltzmann constant, N A represents the Avogadro constant, T indicates the absolute temperature and T ref is the reference temperature. The parameter α a separates the double-well free energy responsible for the phase separation at α a < −4.0 to the homogeneous species distribution without phase separation [22,31]. The parameter α a regulates relative heights of the two wells of free energy. For instance, the relation α a = −2α a means two wells of the free energy at the same energy level. The gradient energy density is expressed as, where κ a is the gradient energy coefficient of the phase a and | • | denotes the norm of a vector. The species diffusion follows the conservation of mass of the form where J denotes the local lithium flux and ∇ · (•) is the divergence of the vector. The species flux is governed by the gradients of the chemical potential, expressed as where μ is the chemical potential and M is the species mobility. The mobility is defined as a function of concentration c of the form [29], where D a is the diffusion coefficient in the phase a. The chemical potential is determined from equation (1) as a functional derivative of the free energy form Note that substituting equations (6)- (8) in the mass conservation equation (5) yields the classical Cahn-Hilliard equation of the form This partial differential equation determines the evolution of concentration c for given initial and boundary conditions. As illustrated in figure 2, a constant concentration in the form of Dirichlet condition is prescribed on one side of the boundary, while on the opposite end a Neumann condition is considered. The remaining two boundaries for 2D (four boundaries for 3D) simulations are considered periodically repetitive. The complete expression can be written in the form, Here c ps is the prescribed constant concentration at the boundary and N x , N y , and N z are the number of grid points in x-, y-and z-directions (x, y, z ∈ x) respectively. Furthermore, a formation of thin layer at the electrode-electrolytes interface is often referred as the electric double layer (EDL). Generally, the thickness of this layer is in the order of only a few nanometers, while active particles are of hundreds to thousands of nanometers scale. This is well reflected in the classical descriptions of the EDL provided by Helmholtz [47], Gouy [48] and Chapman [49]. Preliminary works incorporating EDL with PFM is performed by Guyer et al [50], which was 1D simulations. Despite the intensive research on the EDL recently, majority of works are still limited to 1D due to at least thousands of simulation cells are required to resolve 2-3 orders of length-scale difference between the interfaces and host particles. Typically, the de-/intercalation in the active host material of lithium in Li-ion batteries is approximated by Butler-Volmer type kinetic expressions [51]. Employments of this formulation are still limited to either isolated particles [52] or Newman-type pseudo two-dimensional model (P2D) [53] due to computational complexity and limitations on number of grid points can be selected. Therefore, it is unfeasible to capture such electrode-electrolyte interface effects reasonably even within very powerful techniques like smoothed-boundary method [54,55]. Hence, in the present work, these interface effects are neglected. Instead, an emphasis is placed on differences in diffusional properties of active material and electrolyte. Specifically, the change of concentration dynamics is observed due to differences in diffusities, gradient energy coefficients, and regular solution parameters. In the present study, a standard message passing interface (MPI) based finite-difference scheme is employed on a squared 2D (or cubed for 3D) grid to solve the evolution equation (9) along with the boundary conditions (10). For computational convenience, the variables in the evolution equation and the boundary conditions are normalized as, where L is the reference length (a unit size of the pixel or voxel of the simulation domain), D 1 is the diffusivity (in units of m 2 s −1 ) of lithium species in the electrode particles, and the notation• represents normalized quantities. The spatial derivatives are discretized by the finite difference scheme implemented on a staggered grid, in which the scalar variables are stored in the cell centers, whereas the vectors are located at the cell faces. For instance, the scalar value ∇ · ∇ĉ at discrete time n and spatial location (i, j, k) can be evaluated as, where Δx, Δy, and Δz are the simulation cell widths in x, y, and z-directions, and O(•) represents order of approximation. Besides, the gradients of concentration ∇ĉ in all directions are and the temporal derivatives is discretized by an explicit Euler scheme of the form where Δt is difference between current time step n and next time step n + 1. The simulation grid size employed for 2D microstructures is 200 × 200 cells with approximately seven cells for the interfaces. The number of grid points in a particle is 380, 531, 707, and 1018 pixels for particle sizes R = 0.61, 0.72, 0.83, and 1.00 μm respectively. The wallclock time for a typical 2D simulation to generate microstructure (model described in appendix A) and to study concentration evolution (models described in section 2) are approximately 48 h for each on AMD Ryzen Threadripper 2950X Processor with MPI paralleled code utilizing 16-Cores. While, 3D microstructure of 83 × 83 × 83 simulation cells consumed nearly 103 h of wall-clock time for each. Note that, the choice of boundary conditions does have a significant impact on the wallclock time required for the simulations.

Results and discussion
Simulations are performed in three independent stages. Firstly, the characterisation of electrode microstructure is obtained by Allen-Cahn equation based numerical model described in appendix A. Thereafter, the tortuosity of the obtained microstructure is determined according to the description in appendix B. Finally, Cahn-Hilliard equation based numerical model derived in section 2 is employed to investigate microstructural properties influencing species diffusion in electrode microstructure.
The generic model of species diffusion for multiphase systems of LMO material is described in the above section. As a demonstration of the model, the results are obtained for considering only two phases, the electrode and the electrolyte with highlighting the broad applicability. For instance, the simulation results are presented for propylene carbonate +1 M LiClO 4 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 1.
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 flux (constant concentration condition at the separator). In the present work, we intend to study the discharging process. Therefore, the electrode region is filled with lower equilibrium valueĉ(x,t = 0) = 0.12, while the electrolyte region with its respective equilibriumĉ(x,t = 0) = 0.50. A higher concentration value is prescribed at the boundaryĉ ps (x =N x ,ŷ,ẑ,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 = 100/V V (ψ 1ĉ )dV) frequently. Once there are no changes observed for a reasonable time then the simulations can be stopped. Figure 3 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 Li-rich and Li-poor 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 [36]. Furthermore, it can be depicted from (a3.1) and (b3.1) that the newly exposed interfaces observe higher rate of change of concentration (dĉ/dt) compared to the older ones. 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 phaseseparated traveling front from the separator to the current collector is apparent in the latter case in (b4) for the porous electrode.

Basic features of the model
The porous structure of an electrode is of interest in battery applications due to high surface area yet minimum energy loss [57]. Another type of structure which is widely under scrutiny is the planar electrode, as utilized in solid-state thin film batteries [12]. 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 4. 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 4(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 figures 4(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 3(b1). As the phase separation in a particle starts from the section toward the separator, several particles at the same distance in an affected layer initiate the phase separation simultaneously. Due to enhanced diffusivity of the electrolyte, the lithium flux traverses through the surrounding Here τ denotes tortuosity, ρ indicates porosity of the system, and 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 i ∝ t −1/2 (Cottrell [58] line, where i is the flux of lithium species measured as the rate of change of SOC) 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 JMA equation (16).
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 3(b3). Therefore, the wave-like pattern can be observed in figure 4(b).
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 5(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-and surface-transport limited diffusion [59], which are the standard techniques employed to measure the change in species current in a controlled potential environment in electrochemistry.
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 6(a). Therefore, the transportation rate (or species flux) is controlled by the bulk diffusion. The total time of charge for bulk diffusion limited process observing phase separation behavior is given by a specific time t s . For initial small duration t t s , a simplistic relation of species flux with time in a bulk-transport limited electrode under a potentiostatic condition can be described by Cottrell equation [58], where i (in h −1 units) is the flux of lithium species measured as the rate of change of SOC and k (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) to measure the diffusivity [60]. Therefore, the Cottrell equation (15) for k = 0.75 h −1/2 is plotted with the results obtained from planar electrode for comparison. After the deposition at the particle surface, species transport through the bulk of the electrode. Due to 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 (15), as shown in figure 5(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 = 1.00 μm) tends toward 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 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 6(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 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 [61] relates the transformed fraction in two-dimensional systems as, where ζ denotes transformed fraction and N c is the constant associated with nucleation rate. The JMA equation defines the transformation of the particles to follow a sigmoidal shape. Figure 5(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 reference [59] 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 size of the particles.

Effect of tortuosity
The complete process of phase separation is displayed in figure 3 for two representative systems: (a) electrode particles aligned perpendicular and (b) parallel to the separator. A characteristic time t 80 is defined to investigate transportation dynamics, which signifies time required to observe the 80% SOC. System (a) consumes t 80 = 5.15 h, while system (b) requires t 80 = 13.78 h. As the systems correspond to equal particle sizes and porosity, the distribution of particles can be considered immensely influential on the characteristic time t 80 . One of the parameters that quantify the geometrical particle distribution is the tortuosity of the system τ . Based on studies of references [43][44][45], the technique to measure the tortuosity is described in appendix B. For the provided representative cases in figure 3, tortuosity τ = 1.28 attained t 80 = 5.15 h, while τ = 2.43 necessitates t 80 = 13.78 h. Furthermore, for different τ and maintaining the porosity and the mono-disperse particle sizes of the systems consistent, t 80 is linearly related to tortuosity τ , which is plotted as the red curve in figure 7. Therefore, the increase in tortuosity of the system linearly increases the time required to attain 80% SOC, t 80 .

Effect of particle size
The rate of change of SOC and fraction transformed relate to the particle sizes R, as displayed in figures 5(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 toward 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 bulktransport 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 7 shows the effect of different mono-disperse particle sizes as a function of tortuosity τ and t 80 with maintaining porosity. For different particle sizes, as shown in figure 5, the system of smaller particles reaches saturation earlier compared to larger particles. Apart from that, linearity in the τ and t 80 relation is observed for all systems of various particle sizes. However, τ vs t 80 lines are not parallel for different particle sizes with the same porosity, instead the slope increases with decreasing particle sizes. In addition, the difference in t 80 is higher at lower tortuosity τ compared to the higher ones. Furthermore, the porosity certainly affects the time t 80 , which is described in the next section.

Effect of porosity
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 8 shows the relation of τ and t 80 with the porosity ρ for a consistent particle size R. A linear relation of t 80 with τ can be observed for respective values of porosity.
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 80 . Also, the slope of the linear relationship (τ with t 80 ) is not constant, instead, similar to the observations at different particle sizes, the slope increases with decreasing porosity. Contrary to that observation, the increment in t 80 with change in porosity ρ is lower at lower values of tortuosity τ compared to the higher ones. This implies that the system of higher τ exhibits appreciable increment in t 80 compared to the lower ones for the same porosity difference.
There are some experimental studies which can corroborate the results observed in the presented numerical study and trends reported therein can be compared. For instance, Chen et al [62] reported that the tortuosity of electrode with particles aligned to the lithium transport pathways is lower (≈1.25) than that of perpendicularly aligned electrode (≈4.46), while randomly aligned microstructure shows intermediate value, which is in complete agreement with the present study. Furthermore, high electrode tortuosity is reported to cause thick Li deposition layer [62,63], similar to the layer-by-layer intercalation observed in figure 3(b) of the manuscript. On the other hand, low electrode tortuosity enables homogeneous Li transport and uniform Li deposition across the host [62,63], similar to the concurrent intercalation observed in the present work as shown in figure 3(a). This findings are further supported by reference [64] that the low tortuosity architecture can boost the transport capabilities significantly.
The influence of particle sizes in electrode on lithiation dynamics are reported in references [65,66]. Smaller particle sizes in microstructure maximizing the surface area with electrolytic solution, thus improving transport properties. Besides that, the increase in pore space (porosity) of electrode is reported to provide better specific power and improved lithiation/delithiation in reference [67], which is in agreement with the presented results.

Lithium transport in three-dimensional particle
Finally, figure 9 demonstrates the capability of the employed model that the method can be readily extended to simulate three-dimensional systems. A 3D electrode is considered for the numerical experiment to study the transport mechanism. The potentiostatic condition initiates the formation of the Li-rich phase from the outer surface of the particles, otherwise consists of Li-poor phase throughout, as shown in (b). The Li-rich phase continues to grow further at the surface of the particles in (c)-(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 of phase segregation 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 (τ = 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 our preliminary findings, the depth of the electrode certainly influences the SOC evolution, which will be subject of future investigations.

Concluding remarks
The porous structures are characterized numerically based on Allen-Cahn-type PFM, thereafter phase separation dynamics is investigated with help of Cahn-Hilliard-type formulation on different electrode microstructures. This numerical setup is utilized to understand the effect of various microstructural descriptors such as particle size, porosity, and tortuosity on the lithium transport mechanism in porous electrodes. Our results suggest that the transportation rate of the system is strongly related to the tortuous pathways formed by the particle orientation and distribution, which can be quantified by the microstructure descriptors.
When controlling the lithium concentration at the separator, the lithium transportation rate is observed to be linearly related to the tortuosity. Furthermore, these lines are converging with variation in particle size at higher tortuosity values, while they are diverging with variation in porosity. Furthermore, the obtained numerical results are compared with the bulk-transport and surface-reaction limited theories. The results show that the smaller mono-disperse particles tend to the surface-reaction limited theory, while the larger ones tend toward the bulk-transport limited theory of planar electrode.
A number of interesting directions can be pursued hereafter. For instance, the focus of the present study is on the effect of various geometrical parameters such as monodisperse particle size, porosity, and tortuosity on the charge dynamics. A straightforward extension is to investigate other crucial parameters such as electrode overpotential (value of c ps ), electrode thickness, size polydispersity, and shape of the particles. These parameters have physical significance, which will be considered in future investigations.

Appendix A. Preparation of electrode microstructure
The electrode structure is characterized numerically by employing volume preserved technique on multiphase Allen-Cahn formulation. Each particle is considered as a separate phase freely evolving in the liquid electrolyte phase. A brief procedure is outlined in the following paragraphs.
The inhomogeneous system consists of different phases whose identity is distinguished by its physical state or orientation. In the present description, a multi-phase Allen Cahn model is considered to study the evolution of multiphase systems whose interfacial energy decreases with preserving their volume [68]. The evolution of a general system containing N phases is governed by Ginzburg-Landau free energy, where f AC (φ, ∇φ) is the free energy density, V is the domain under consideration, a vector phase-field parameter φ( x represents spatial coordinates and denotes the thickness of the diffuse interface. The gradient energy density is expressed as [69], where γ αβ denotes the interfacial energy density between α and β phases, q αβ = φ α ∇φ β − φ β ∇φ α , which defines the gradient vector in the normal direction to α − β interface. The function a αβ (φ, ∇φ) represents the interface anisotropy. For the present study, a faceted anisotropy is considered of the form, where {η i,αβ |i = 1, . . . , η αβ } represents position vectors in the Wulff shape of phase α concerning phase β. The bulk energy density functional w AC (φ) is assumed to be a multi-obstacle type potential, The higher-order interfacial energy term γ αβδ penalizes the ghost phase occurrences at the interfaces between two phases. The additional bulk free energy density g AC (φ) is responsible for the volume preservation, can be expressed as, where the function h(φ α ) = φ 3 α (6φ 2 α − 15φ α + 10) interpolates the free energy density terms χ α between the bulk phases. The time-dependent antiforce terms χ(t) ∈ {χ 1 , . . . , χ N } are calculated at each timestep to maintain the phase volume equal to initially prescribed. For a detailed algorithm, the readers are referred to reference [70]. Finally, the evolution of the phases is governed by a general model for multi-phase Allen-Cahn type equation as where τ R denotes the relaxation coefficient, ∂(•)/∂t is partial derivative with time t and δ(•)/δφ α denotes functional derivative. The second term in the above equation ensures the constraint N α=1 φ α (x, t) = 1. Inherently, the isolate boundary condition implies ∇φ α · n = 0, where n is the normal vector at the wall of the simulation boundary. In other words, the phases form a 90 • angle at the boundary. To specify other contact angles, an extension is made in the form of ) denotes the boundary function and the terms γ αb represent the interfacial energy density between the α phase and the boundary, the value of these parameters controls the contact angle between the boundary and different phases from 0 • (complete non-wetting) to 180 • (complete wetting). The term in the parentheses in equation (A.7) ensures the sum of all phases at any location in the simulation domain is unity.
The electrode microstructure is obtained from the random distribution of N − 1 circular (or spherical in 3D) particles in the domain of electrolyte. Each particle corresponds to a separate phase φ α and the last phase φ N denotes the electrolyte. These phases are allowed to evolve under N equation (A.6) with volume preserved for non-wetting boundary conditions (A.7). The particles evolve to a geometry specified by the position vectors η i,αβ to the vertices in equation

Appendix B. Calculation of tortuosity
The electrode particles are surrounded by the liquid electrolyte, which diffuses lithium species in the electrode. The species diffusivity of the electrode is much lower than the electrolyte [56]. Also, the insertion and the extraction of species in the electrode particles takes place at the interface between the electrode and the electrolyte. Therefore, understanding of tortuous pathways in the electrolyte is of scientific interest. In the present study, instead of assuming a Bruggeman relation [42], a computational tool is employed to calculate the tortuosity of the electrolyte microstructure considering infinite resistive electrode particles [43].
The tortuosity is calculated in the vertical direction, where the species flux propagates, from the separator at the top (Φ = 1.0) to the current collector at the bottom (Φ = 0.0). A potential difference Φ ∞ is applied between two ends, which is predefined as 1.0 V at the separator (top) and 0.0 V at the current collector (bottom) for the present investigation as shown in figure A.10(c). While periodic boundary conditions applied on the remaining two (four for 3D) boundaries. The potential distribution inside the simulation domain is calculated using the Laplace equation of the form, where σ ety is the electrical conductivity of the electrolyte. Afterward, the current density i = σ ety ∇Φ can be estimated at each pixel (voxel in 3D). The total current flow from the cross-section perpendicular to the x-direction can be obtained as a summation, where A is the area of the cross-section surface under consideration. Note that, due to conservation, the measured current flow at any cross-section in horizontal direction should be nearly equal. If not then the meaningful conclusions can not be drawn from the presented method. For instance, in the case where there is no path for the species to flow across when electrode particles completely block the path by forming an entire blockage in the horizontal direction.
The tortuosity is a measure of conductive pathways deviated from the ideal straight channels of uniform cross-section. Therefore, the tortuosity estimation is based on the actual current flow in the microstructure in contemplation to the current flow in the ideal path without any discontinuities, where I ideal is the ideal current flow, determined from the analytical expression I ideal = σ ety Φ ∞ A/l and l is the length of the domain in the vertical direction.

Appendix C. Phase separation in electrode along with homogeneous mixture in electrolyte
The behavior of two-phase coexistence in the electrode compounds are well known during the insertion and the extraction of lithium species [13][14][15]17]. However, instead of phase separation, a homogeneous distribution is apparent in the electrolyte. Both the behaviors can be achieved by careful consideration of the free energy parameter α a in equation (3) and the gradient energy parameter κ a in equation (4). The Cahn-Hilliard equation is employed for the study of phase separating behavior in ample literature [22,[25][26][27][28][29][30]. Therefore, we only demonstrate that for a specific choice of these parameters, the Cahn-Hilliard equation exactly corresponds to the Fick's dilute solution model [71]. Considering a special case, where any phase a consists of the regular solution parameter α a = 0 and the gradient energy parameter κ a = 0. The chemical potential of such a system can be expressed as, Here μ a is the chemical potential in the phase a. The species flux can be written as, Here J a and M a represent the species flux and the mobility in the phase a respectively. Equation (C.2) is equivalent to the equation responsible for the homogeneous mixture. Therefore, parameters α a = 0 and κ a = 0 can be utilized for the regions exhibiting non-phase separating behaviors, such as in the electrolyte. This can be identified in figure C.11, where phase separation is observed in electrode particles, while non-phase separating behavior is depicted for electrolyte.