Microstructure-informed analytical modeling of composite cathode for solid-state batteries

Experimental LCO-LLZTO composite preparation & imaging

LiCoO₂ (MTI Corporation) and Li_6.4La₃Zr_1.4Ta_0.6O₁₂ (NEI Corporation) powders were ball-milled together in a ratio of 75 wt%/25 wt% using a SPEX 8000D Mixer/Mill. The powder mixture was then pressed into a pellet at 3000 psi for 2 min using a ¹/₂” die (MTI Corporation) and a hydraulic press (YLJ-15, MTI Corporation). A small amount of acetone (40 μL) was added to the powder mixture inside the die before pressing to help with pellet release. The pellet was then sintered in a zirconia tray in a box furnace (MTI Corporation) at 1,050 °C in air for 2 h (with a heating and cooling ramp rate of 5 °C/min). It was then broken with tweezers, and the cross-section was analyzed by scanning electron microscopy (SEM) (Thermo Fisher Scientific Apreo).

Simulation workflow

The simulation workflow consists of four steps, as illustrated in Figure 1. First, three-dimensional virtual microstructures were generated using a custom-built code, StructGen, which efficiently constructs multiphase porous microstructures with controlled phase population, particle sizes, and porosities by stochastically populating spherical particles within a simulation domain with periodic boundaries. Overlapping between particles of the same phase or between different phases can be made permissible or impermissible. In the present study, we generated two-phase porous microstructures consisting of two solid phases (Phases 1 and 2). Second, geometric measurements were performed on the generated virtual microstructures to estimate the average domain size of each solid phase ($$ \underline d_{CAM} $$ and $$ \underline d_{SE} $$) and the specific interface area between them (A_s). Third, flux-based simulations were conducted by numerically solving the steady-state diffusion equation under the assumption of Fickian diffusion to compute the effective lithium-ion diffusivity (D_eff) and determine the tortuosity factor (τ²). Finally, the calculated D_eff, τ², and A_s as functions of phase fractions and porosity were analyzed and fitted using analytical models based on a GEMT model for D_eff and τ² and the Dirichlet distribution for A_s. The percolation thresholds in terms of the density (ρ_c) and the SE phase fraction (ρ_c’) were determined from the analytical model.

Figure 1. Overview of the workflow for developing surrogate models to evaluate the effects of phase composition on lithium-ion transport in composite cathodes through microstructure-aware modeling. SE: Solid electrolyte; CAM: cathode active material.

Microstructure generation

The simulation box was discretized into a voxel-based grid of size 200 × 200 × 200, corresponding to a physical size of 20 μm × 20 μm × 20 μm, with each voxel representing a 100 nm cube. The virtual microstructures were created in three steps, as illustrated in Figure 2. For visualization clarity, only 2D cross sections of the 3D microstructure were shown.

Figure 2. Workflow for generating virtual microstructures. A 2D cross-section of the 3D microstructure is displayed for enhanced visualization clarity. SE: Solid electrolyte; CAM: cathode active material.

First, a porous structure of Phase 1 was constructed by randomly populating spherical particles within the simulation box until a targeted phase fraction ϵ₁ was achieved. The particles were modeled as spheres with diameters sampled from a Gaussian distribution (mean: 30 voxels or 3 μm, standard deviation: 5 voxels or 0.5 μm). Overlapping between Phase 1 particles was permitted to ensure efficient packing. Second, the pore regions were filled by inserting Phase 2 particles, also modeled as spheres. Instead of randomly populating Phase 2 particles all at once, an iterative algorithm was proposed to more efficiently achieve the targeted porosity ϵ₃. Specifically, a new Phase 2 particle was inserted during each iteration at the location within the pore phase where the largest distance from the Phase 1 surface was identified based on Euclidean distance search. The maximal distance was used as the radius of the newly inserted Phase 2 particle, ensuring interfacial contact between Phase 1 and Phase 2. This approach essentially requires that each Phase 2 particle corresponds to the maximal inscribed sphere of the pore phase at each iteration. Using this method, the targeted porosity ϵ₃ was quickly achieved within a few iterations, while overlaps between Phase 1 and Phase 2 particles were strictly prohibited. The resulting volume fraction of Phase 2 is determined to be ϵ₂ = 1 - ϵ₃ - ϵ₁. Third, a geometric sintering process was applied using an image filter to ensure percolation and connectivity within the Phase 2 network. This process fused neighboring Phase 2 particles by replacing each voxel of the pore phase with its most common neighboring voxel type. The geometric sintering step results in a slight reduction (approximately 1%) in actual porosity compared to the targeted value but ensures the connectivity of the Phase 2 network.

Following this procedure, a total number of 258 distinct 3D microstructures of the two-phase porous composite were generated, with varying phase fractions and porosities. The resulting ranges of phase volume fractions were as follows: 10%-85% for CAM, 1%-89% for SE, and 0.3%-87% for pores. The statistical distributions of the total phase fractions (ϵ_CAM, ϵ_SE), phase fractions with respect to the solid phase (ϵ_s,CAM, ϵ_s,SE), porosities (ϵ_pore), and density (ρ = 1 - ϵ_pore) are provided in Supplementary Figure 1. The two types of phase fractions are related by ϵ_CAM = ρ × ϵ_s,CAM and ϵ_SE = ρ × ϵ_s,SE.

Geometric measurement

As overlapping between Phase 1 particles was permitted during the microstructure generation, direct segmentation to determine the particle size distributions is not feasible. Instead, the domain size of the CAM and SE phases ($$ \underline d_{CAM} $$ and $$ \underline d_{SE} $$) was used as a proxy to represent the average size of each phase. For a two-phase structure, the domain size is $$ \underline d $$ defined as the reciprocal of the first moment of the structural factor. This quantity can be conveniently calculated using the Fourier transform of the shape function corresponding to each phase, following the method described in Ref.29^[29].

The specific interface area between the two solid phases was calculated as the difference between the interface areas of each solid phase with respect to the pore phase, expressed as: $$ A_{s}=A_{s}^{S E}+A_{S}^{C A M}-A_{S}^{ {pore }} $$. Here, a mask function was applied to isolate each phase, and the ComponentMeasurement function was used to compute the location and number of voxels at the surfaces. All algorithms for microstructure generation and geometric measurements were implemented in Wolfram Mathematica.

Effective property calculation

The effective property calculation has been widely used to characterize transport properties and evaluate the tortuosity factor of complex microstructures in battery components, such as porous separators and electrodes of lithium-ion batteries^[31,32], polycrystalline^[30], or porous solid electrolytes^[33], and composite cathodes for ASSB^[20,28]. This approach is commonly referred to as the flux-based simulation, the Fickian diffusion model, or the direct numerical simulations^[20,31]. It numerically solves the Poisson equation with a spatially varied transport coefficient under an applied gradient of a scalar field variable along each spatial direction. The effective transport coefficient is derived from the resulting average flux density over the applied gradient. In terms of the Fickian diffusion model, the governing equation can be written as

where c(x) is the concentration of the transport species, D(x) is the spatially dependent diffusivity tensor, and x is the position vector in 3D space, i.e., x = (x₁, x₂, x₃). An applied boundary condition imposes a fixed concentration gradient ∇c along one of the spatial directions. The effective diffusivity tensor, D_eff, is then determined as

where <•> stands for spatial average over the simulation box, and the flux density J(x) is calculated as

While the local diffusivities of CAM and SE phases are anisotropic in general, we adopted an isotropic approximation for simplicity. The resulting effective diffusivity D_eff is also a second-rank tensor, and we averaged the diagonal components to calculate D_eff, i.e.,

This approximation is justified in that the virtual microstructures were generated using a stochastic approach, and there is no strong directional preference for the conductive networks. The tortuosity factor, τ², is calculated using the relationship

where D_SE is the diffusivity of the conductive SE phase, D_eff is the effective diffusivity of the composite, and ϵ_sol is the volume fraction of the solid phase, i.e., the CAM + SE phases in the present study. This field-based approach to estimating τ² has shown good agreement with the microscopic method using random walk models and Monte Carlo simulations^[23,34-36].

The governing Equation (1), along with the boundary conditions, is numerically solved using the Fourier-spectral iterative perturbation method^[37]. Specifically, we employed a highly efficient, custom-built FFT-based numerical solver, MesoMicro, which has been benchmarked and detailed in prior studies^[29,30].

Parameterization of diffusivities

To calculate the effective diffusivity of the composite, the diffusivities of SE and CAM phases were parameterized based on values reported in the literature. In this study, we considered the diffusivities of composite cathode active materials LCO, NMC111, NMC622, NMC811, and lithium iron phosphate (LFP) and SE materials with diffusivities similar to LLZO, such as the thiophosphate family. Each of these materials is significant in all-solid-state battery research, and their transport properties have been reported in prior studies. We emphasize that the modeling framework is general and can be easily applied to other combinations of SE and CAM materials, as the results depend only on the ratio of diffusivities in the CAM and SE phases.

In general, the SE has a higher lithium-ion diffusivity than the CAM. The literature, based on theoretical calculations and experimental measurements, indicates a wide range of diffusivities for the active materials^[38] and even within a single SE, such as LLZO. The diffusivity within a single CAM depends on lithiation states; for example, LCO’s diffusivity can range over four orders of magnitude depending on lithiation^[39-41]. The wide range of diffusivities in LLZO depends on differences in processing conditions and physicochemical factors such as local atomic configurations^[42], dopants^[43], interfaces^[30,44], and lithium concentrations^[45]. The diffusivity of the SE was set to 10^-8 cm²/s, based on the atomic simulation-informed study that accounted for the microstructural impacts of polycrystalline LLZO^[30,42]. Many families of solid electrolytes, such as thiophosphates and perovskites, exhibit diffusivities between 10^-11 and 10^-8 cm²/s and approach diffusivities of 10^-7 m²/s for the most optimized materials^[39-41]. For this study, the diffusivity of the CAM ranged from 10^-12-10^-9 cm²/s, which covers the diffusivity of widely researched NMC111, NMC622, NMC811, LCO and even LFP^{[38,40,46,47]}. The highest diffusivity of 10^-9 cm²/s is based on a three-dimensional phase-field study that accounts for lithium intercalation-induced stress effect in polycrystalline LCO^[48]. The diffusivity of the pore phase was set to zero. Since only the relative ratio of diffusivities between the two materials was required, the diffusivity of the SE was normalized to unity. The maximum effective diffusivity corresponds to the limit of a fully dense SE. To account for the variations of CAM diffusivity, sensitivity tests were performed using ratios of D_CAM/D_SE = 0.1, D_CAM/D_SE = 0.01 and D_CAM/D_SE = 0.001. Considering that our effective diffusivity model assumes linear diffusion laws, these ratios represent various combinations of specific CAMs and SEs. Only the largest ratio is presented in many of the main document figures, while ternary diagrams for the two smaller ratios are given in Supplementary Figure 2 and will be analyzed in the Discussion section. We discuss how the percolation threshold depends on this ratio.

Generalized effective medium theory

The effective medium theory (EMT) is a widely used model for predicting the transport properties of heterogeneous materials, particularly for two-phase composites. It provides an analytical approximation for the effective transport properties of a two-phase composite by accounting for the local properties and volume fractions of the individual phases. This approach has been extended to multiphase systems^[33,49], known as GEMT.

We proposed a GEMT model for the ternary system of porous CAM-SE by combining two EMT models of two binary subsystems. The first EMT model describes the interaction between Phase 1 and Phase 2, forming a lumped solid composite phase. The second EMT model then accounts for the mixing between the pore phase and the lumped solid phase. This nested approach enables the calculation of the overall effective diffusivity, D_eff, through the formula

where ρ_c is the percolation threshold of the three-phase composite in terms of the solid phase fraction (i.e., density), and t is the effective percolation slope^[33,49,50]. The lumped solid phase has a diffusivity, D_solid, which is calculated based on the diffusivities and volume fractions of the individual SE and CAM phases, i.e.,

where ρ_c’ is the percolation threshold between SE and CAM, and t’ is the corresponding effective percolation slope. By solving Equations (6) and (7), one can obtain an analytical expression of D_eff in terms of ρ and the solid volume fraction ϵ_s,SE, i.e., D_eff(ρ, ϵ_s,SE), with four fitting parameters t, t’, ρ_c, and ρ_c’. To be consistent with our numerical simulation, normalized values of D_SE and D_CAM are used and set to be 1.0 and 0.1, respectively, unless otherwise stated. A negligibly small value (10^-7~10^-4) is assumed for the normalized D_pore to aid in data fitting, depending on the choice of D_CAM/D_SE.

Representative microstructures

We fabricated a simple LiCoO₂-Li_6.4La₃Zr_1.4Ta_0.6O₁₂ (LCO-LLZTO) composite by ball milling 25 wt% LLZTO and 75 wt% LCO powders together for 60 min, pressing the mixture into a pellet, and co-sintering it at 1,050 °C in air for 2 h. We used SEM to characterize the microstructures of the cross-sectional samples, as shown in Figure 3A. The backscattered SEM reveals a relatively dense structure with a clear distribution of two phases: LCO, represented by dark gray regions, and LLZTO, represented by light gray regions. Microporous regions, represented by the dark voids, are observed near some interfacial areas, while smaller voids are also seen inside the bulk phases. Both phases are composed of interconnected particles that are well-sintered, forming intertwined pathways for Li-ion and electronic transport. To quantify the phase populations and porosity, we segmented the two phases based on their grayscale contrast within the selected rectangular region outlined in dashed white in Figure 3A. The analysis determined volume fractions of approximately 61% LCO, 28% LLZTO, and 11% porosity, which agrees reasonably well with the initial mass fraction of the powders (25 wt% LLZTO and 75 wt% LCO). Note that at 61% CAM, it is likely to form an electronic percolating network and thus electronic conduction, provided that the CAM is also a good electronic conductor. In reality, however, a conductive carbon additive and binder are often added to composite cathodes to improve electronic conductivity^[51]. For the LCO-LLZTO composite, co-sintering can eliminate the need for a binder, while a conductive agent would be added to increase the electronic conductivity, which is not explicitly included here. Instead, we focused on a simplified system to explore the effects of the phase fractions of these two main components (LCO and LLZTO).

Figure 3. (A) Backscattered SEM cross-section of a co-sintered 75 wt%/25 wt% LCO-LLZTO composite. Inset: the phase fraction is estimated to be around 61% LCO, 28% LLZTO, and 11% porosity based on the grayscale histogram of the selected rectangular region outlined in dashed white. The black, dark gray, and light gray regions correspond to voids, LCO and LLZTO, respectively; (B) A representative example of microstructure generated stochastically from spherical particles. The transparent, red, and yellow domains correspond to pores, CAM, and SE phases, respectively. SEM: Scanning electron microscopy; LCO-LLZTO: LiCoO₂-Li_6.4La₃Zr_1.4Ta_0.6O₁₂; CAM: cathode active material; SE: solid electrolyte.

Inspired by the experimental microstructure, we constructed 3D virtual microstructures to qualitatively represent and capture the composite cathode microstructure using the stochastic approach described in Figure 2. An example microstructure is shown in Figure 3B wherein the CAM phase is colored in red, the SE phase in yellow, and the pores are left transparent. Pores near the SE-CAM interface and CAM grain boundaries are observed, consistent with the experimental observation. The volume fractions of the three phases (64.3% CAM, 23.1% SE, and 12.5% porosity) compare well with the experimental measurements. However, the pores in the 3D microstructure appear more dispersed compared to the localized voids seen in the SEM cross-sections. In general, the virtual microstructure successfully captures the particulate morphology of the composite microstructure observed in the SEM image.

Instead of a direct, one-to-one comparison between the experimental and simulated microstructures, our focus is on exploring how variations in phase fractions and porosity influence the geometric and physical properties of the 3D microstructures. Therefore, a few representative microstructures with different phase fractions and porosities are shown in Figure 4. Figure 4A represents a highly porous microstructure with comparable fractions of CAM and SE phases. The particle size of the SE phase is also set to be larger than that of the CAM phase, resulting in a dispersed distribution of CAM particles, which may restrict the electronic conduction. Figure 4B shows a more ideal dense structure with a dominant CAM phase and comparable particle sizes. The resulting D_eff and A_s are expected to be higher than the structure in Figure 4A. With a slight increase in targeted porosity compared to Figure 4B, Figure 4C illustrates a more realistic structure with a moderate amount of pore phase while maintaining an approximate 3:1 ratio between CAM and SE phases. It is expected that the D_eff and A_s in Figure 4C should fall between those in Figure 4A and B.

Figure 4. Representative 3D microstructures within the 250 dataset with phase fractions labeled. (A) A nonideal microstructure where the SE phase fraction is too low and the particles of both phases are isolated; (B) An ideal highly dense microstructure close to the predicted SE percolation threshold; (C) A microstructure with phase fractions and porosity closest to the experiments. SE: Solid electrolyte; CAM: cathode active material.

Numerical simulation

To understand the impact of phase fractions on the geometric and physical properties of the composite microstructures, we present the results of effective diffusivity (D_eff), reciprocal of the tortuosity factors (1/τ²), and specific interface area (A_s) as functions of the phase fractions of CAM, SE, and pores, represented by ternary diagrams in Figure 5. The original data are plotted as dots without interpolation. Each data point is colored according to its hue value, as indicated by the associated color scheme. In total, 258 data points are shown for each plot.

Figure 5. Ternary diagrams of (A) effective diffusivity (D_eff), (B) reciprocal of tortuosity factor (1/τ²), and (C) specific interface area (A_s). The total volume fraction of solid electrolyte (ϵ_SE) is on the horizontal axis (in red), the fraction of cathode active material (ϵ_CAM) is on the vertical left-hand axis (in green), and the porosity (ϵ_pore) is on the vertical right-hand axis (in black). The labels (a), (b), and (c) in (A) correspond to the microstructures in Figure 4A-C, respectively. There are three auxiliary lines to aid in reading the phase fractions of Point (c) in (A), where the dashed green line indicates 64% CAM, the dashed black line indicates the 13% porosity, and the dashed red line indicates 23% SE. SE: Solid electrolyte; CAM: cathode active material.

Figure 5A shows that D_eff is highest near the SE apex due to the higher diffusivity of the SE versus the CAM. Reducing the SE fraction or increasing the porosity leads to a decrease in D_eff, with porosity having a more significant influence. For instance, at 50% porosity, D_eff drops almost to zero, whereas a finite D_eff remains even near the CAM end. For 1/τ², which is derived from D_eff and density (see Equation 5), its trend generally aligns with that of D_eff, albeit with slight differences toward the CAM apex. Specifically, while D_eff drops below 0.2 near the CAM apex, the reciprocal of the τ² decreases more gradually, reaching only half its value at the CAM apex.

In contrast, A_s exhibits a distinct trend in the ternary diagram compared to the other two quantities [Figure 5C]. First, the diagram is almost symmetric with respect to the phase fractions between SE and CAM because A_s is a geometric factor independent of the relative magnitude of diffusivities. Second, the maximum A_s is observed at a 50%:50% intermixing of SE and CAM with zero porosity. If porosity increases slightly, A_s drops sharply. This finding is consistent with previous studies on two-phase mixtures with similar particle sizes^[13].

Comparison of Figure 5A-C reveals clear differences and correlations in how composition affects D_eff, τ², and A_s. Geometric factors of the 3D microstructure (τ² and A_s) provide valuable insights into structural arrangement but cannot fully capture the physics-derived effective property (D_eff). For instance, while D_eff approaches negligibly small values near the CAM apex, the tortuosity factor (τ²) remains at an intermediate value. This behavior is understandable because the CAM is assumed to be a Li-ion conductor in our study. As a result, the CAM phase contributes geometrically to tortuosity, even though its diffusivity is limited. This distinction underscores the need for flux-based simulations rather than relying solely on geometric descriptors such as τ² and A_s, as these simulations directly incorporate the interplay between physical ion transport and the microstructural connectivity.

Moreover, Figure 5 highlights a tradeoff in optimizing the phase compositions of composite cathodes to balance ion transport and interfacial reaction kinetics. While increasing SE volume fraction leads to enhanced mass transport, intermixing of SE and CAM with a similar phase fraction is essential for sufficient interfacial reactions. Therefore, it is of critical importance to integrate both geometric and transport-based considerations in the design process of composite cathodes. In reality, increasing the SE often improves both the mass transport and intermixing, since the phase fraction of SE is almost never above 50% because the fraction of active material is important for energy density.

Analytical models

To construct a surrogate analytical model of D_eff, we employed a modified GEMT tailored for three-phase composites with percolation effects between each two components. This formulation combines two EMT models, each representing one of the binary subsystems, as detailed in the Methods section. By solving Equations (6) and (7), we obtain an analytical expression of D_eff(ρ, ϵ_s,SE) as a function of the total density ρ and the solid-state volume fraction of SE ϵ_s,SE.

The analytical model is governed by four fitting parameters: the two percolation thresholds (ρ_c and ρ’_c) and the percolation slopes (t and t’), corresponding to the respective EMT formulations. These parameters were determined by nonlinear regression using a least-squares objective function to fit the simulated D_eff(ϵ_CAM, ϵ_SE). The conversion between the overall volume fractions and solid-state volume fractions was accounted for in the fitting procedure. This approach enables the estimation of key features, such as percolation thresholds and scaling exponents, which are essential for developing accurate surrogate models of effective transport properties in multiphase composites.

The regression analysis yields fitting parameters ρ_c = 0.38 and t = 1.8 and ρ_c’ = 0.33 and t’ = 0.89. These values define the percolation thresholds between the pore-solid phases and the SE-CAM solid phases, respectively. According to classical percolation theory, the effective conduction coefficient (D_eff in our case) scales as D_eff ∝ ( - ρ_c)^t and D_eff ∝ (ϵ_s,SE - ρ_c’)^t’ above the respective thresholds. The exponents (t and t’) characterize the sharpness of the transition from non-percolating to percolating regimes. Figure 6 compares the analytical model predictions with the simulation data of D_eff. The close agreement confirms that the GEMT-based surrogate model captures the key trends across the entire composition space.

Figure 6. Data fitting using the GEMT model of D_eff as a function of (A) the density ρ and (B) the solid phase volume fraction ϵ_s,SE; (C) A perspective view of the simulation data points and the functional surface by the analytical model over the ternary composition domain. GEMT: Generalized effective medium theory; SE: solid electrolyte.

To visualize the model fit and dissect the contributions of different percolation mechanisms, we projected the simulation results along two axes, as shown in Figure 6A and B. In Figure 6A, D_eff is plotted against the total density ρ. The model predictions accurately bound the simulation data between the theoretical upper and lower limits, corresponding to single-phase SE (ϵ_s,SE = 100%) or CAM (ϵ_s,CAM = 0%), respectively. A sharp transition point near ρ = 0.4 is seen, in good agreement with the fitted predicted threshold ρ_c = 0.38. This result suggests that a minimum density of approximately 62% is required for percolation in the solid matrix.

In Figure 6B, the same dataset is reprojected along the ϵ_s,SE axis to isolate the transport transition between SE and CAM phases. Here, we highlight data corresponding to porosities in the range of 0%-15% and render higher-porosity samples in gray, as such structures are less likely to exhibit continuous transport pathways and are less relevant to dense cathode architectures. In fact, experimental studies on composite cathodes (e.g., LCO-Li₆PS₅Cl) report porosities of 20%-30% in uncalendered form, which can be reduced to 13% through mechanical pressing, and further down to 6%-8% via cold pressing^[14]. The analytical model accurately captures the variation of D_eff within the low porosity regime, including a crossover from CAM-dominated conduction and SE-dominated transport. This transition corresponds to a percolation threshold between the two solid phases, with the fitting suggesting a critical SE content of approximately 33% for solid-solid percolation.

Comparison with theoretical limits from 3D continuum percolation theory suggests that our estimated thresholds are slightly elevated, while the exponents are being underestimated, with classical theory predicting ρ_c = 0.29 and t = 2 for for conductivity^[52,53]. The deviations are attributed to factors such as the finite diffusivity of the CAM phase, permissible particle overlap, finite-size effects of the particle packing models, and numerical artifacts in the calculation of D_eff. Nevertheless, the fitted model effectively captures the essential percolation behavior observed in simulations and provides a robust analytical framework for rapid evaluation of transport properties. Moreover, we found that the fitting parameters vary systematically with the assumed diffusivity ratio between SE and CAM, with the percolation threshold ranging from 0.25 to 0.33 as the contrast varies from 0.001 to 0.1. This sensitivity and its implications will be further discussed in the next section.

The active interfacial area A_s between the CAM and SE is a critical microstructural descriptor that governs the efficiency of lithium-ion intercalation. Maximizing this interfacial area is essential for enhancing reaction kinetics and ensuring uniform current distribution. To model the dependence of A_s on the phase composition, we employed an analytical expression derived from the probability density function of a Dirichlet distribution of order three. The model is expressed as

where Г(·) is the Gamma function, α₁, α₂, and α₃ are shape parameters, and k is a scaling factor. Using nonlinear least-squares regression, we identified the best-fit parameters as α₁ = 2.3, α₂ = 2.2, α₃ = 0.93, k = 0.011.

Figure 7A presents the variation of A_s with respect to the solid-phase volume fraction of SE ϵ_s,SE. The model accurately reproduces the peak in the interfacial area near an equal-volume mixture of CAM and SE, corresponding to optimal CAM-SE contact. Furthermore, Figure 7B visualizes both the simulation data and analytical predictions in a ternary 3D plot, showing strong agreement across the full compositional space.

Figure 7. (A) Scatter plot of calculated specific interface area (A_s) as a function of the solid volume fraction of SE with an upper bound fit using the Dirichlet distribution function (black dashed curve); (B) A perspective view of the simulation data points and the functional surface by the analytical model over the ternary composition domain. SE: Solid electrolyte; CAM: cathode active material.

These findings underscore the critical role of high solid-phase packing density and low porosity in promoting interfacial contact. The presence of voids reduces the effective interfacial area, thereby limiting the number of active reaction sites and compromising electrochemical performance. The Dirichlet-based analytical model developed here provides a computationally efficient and accurate method to estimate A_s across arbitrary phase compositions, offering a valuable tool for microstructural optimization in solid-state battery electrode design.

Figure 8 presents the calculated D_eff, 1/τ², and A_s using the three parametrized analytical models in terms of the overall phase fractions and porosity. These ternary diagrams enable quick identification of corresponding values for any given combination of phase fractions and porosity, aiding in the design of composite cathodes to optimize effective diffusivity while preserving sufficient active interfaces. As an example, three points labeled (a), (b), and (c) are highlighted in Figure 8A, corresponding to the three microstructures in Figure 4A-C, respectively.

Figure 8. Ternary diagrams of (A) the effective diffusivity D_eff, (B) the inverse of tortuosity factor (1/τ²), and (C) the specific interface area (A_s) using corresponding analytical surrogate models. The parameters for the GEMT are t’ = 0.89, t = 1.75, ρ_c = 0.39, ρ_c’ = 0.33 and the parameters for the third-order Dirichlet distribution is α₁ = 2.3, α₂ = 2.2, and α₃ = 0.93. The labels (a), (b), and (c) in (A) correspond to the microstructures in Figure 4A-C, respectively. There are three auxiliary lines to aid in reading the phase fractions of Point (c) in (A), where the dashed green line indicates 23% CAM, the dashed black line indicates 13% porosity, and the dashed red line indicates 64% SE. SE: Solid electrolyte; CAM: cathode active material.

To determine the phase fraction for a specific point on the ternary diagram, auxiliary lines parallel to the three axes can be drawn through the point, and the values can be read directly from the respective axes. For instance, Point (c) in Figure 8A corresponds to 23.1% SE, 64.3% CAM, and 12.5% porosity. A horizontal dashed black line is drawn parallel to the axis ticks of the pore phase, yielding a value of 0.13 from the axis. Similarly, a dashed green line parallel to the CAM axis ticks gives a value of 0.23, while a dashed red line parallel to the SE axis ticks gives a value of 0.64. The corresponding D_eff value can be read off from the color scheme. A similar analysis can be performed for the other points. By comparing the color schemes of the three points, it is evident that D_eff follows the trend Point (b) > Point (c) ≈ Point (a), which aligns with the visual intuition from their microstructures in Figure 4. This consistency confirms the rationality of the analytical modeling results.

Similar analyses can be conducted to evaluate τ² and A_s for a given phase fraction and porosity or to select optimal phase fractions that balance overall diffusivity and interfacial area at a fixed porosity. Comprehensive use of Figure 8A-C provides an efficient and effective tool for designing composite cathodes, enabling consideration of the tradeoff between bulk transport and interfacial reaction kinetics.

Factors influencing the simulation results

In the following, we discuss the differences between our results and those of prior studies regarding assumptions about ionic transport in the composite cathode, the definition of the tortuosity factor, and the general applicability of our results due to the flexibility of our model and sensitivity tests that vary the Li-ion diffusivities in the CAM.

First, in most prior studies, the CAM was typically treated as an electronic conductor and ionic insulator, while the SE was considered an ionic conductor and electronic insulator^[5,6,19,31]. In contrast, our methodology assumes that both CAM and SE are conductive to lithium ions. This assumption is supported by the non-negligible diffusivities of lithium ions in many CAMs at room temperature, which facilitates effective mass transport in both the CAM and SE. This fundamental difference in assumptions has implications for the interpretation of transport properties and the estimation of related parameters such as the tortuosity factor τ².

Second, the estimation of the tortuosity factor (τ²) in our study differs from conventional approaches. In prior work, the tortuosity factor was often calculated with non-percolating dead-end pores removed to enhance numerical stability and computational efficiency. In contrast, our approach does not exclude these components and therefore implicitly accounts for their contributions to D_eff. This approach is partially justified by the suggestion to use the electrode tortuosity factor (τ_e²) instead of the conventional tortuosity factor (τ²) for more accurate measurement of the transport properties of porous electrodes, as dead-end pores have been shown to play a critical role^[54]. Moreover, tortuosity is inherently anisotropic^[32], and in our study, the tortuosity factor can be derived from the anisotropic effective diffusivity tensor D_eff. However, due to the random nature of microstructure generation in our simulations, there is no strong directional preference for ionic conduction, resulting in nearly isotropic tortuosity factors. It is important to note that while tortuosity factors have been extensively studied in the context of porous electrodes and separators in lithium-ion batteries, their definition and significance in composite cathodes for ASSB remain less well understood and deserve further investigation.

Finally, we performed a sensitivity test to investigate how the diffusivity in the CAM, especially compared to the SE, affects the effective diffusivity and percolation thresholds of the overall density and fraction of SE. In our study, we set D_LCO/D_LLZO = 0.1 for data presented in the main text, but also performed sensitivity tests using D_LCO/D_LLZO = 0.01 and 0.001. For each set of data, we refit a corresponding GEMT model, as shown in Supplementary Figure 3, and extract the best-fit parameters, as listed in Table 1.

The fitting parameters show a trend: as D_CAM decreases relative to D_SE, the density needed to reach the percolation threshold increases, as expected. Conversely, the fraction of SE needed to reach the percolation threshold decreases from 33% to 25%, suggesting that a smaller SE fraction is sufficient to shift the diffusion behavior toward SE pathways. This observation is consistent with the underlying physics: when the CAM diffusivity is much lower than that of the SE, ions preferentially use SE percolative pathways once they form. In contrast, if the CAM diffusivity is much closer to that of the SE, the ions will use the larger volume fraction CAM pathways before switching over to the SE percolative pathways as that volume fraction increases. Plotting the data with respect to the solid-phase volume fraction of SE ϵ_s,SE [Supplementary Figure 4] shows that D_eff exhibits different slopes before reaching the percolation threshold ρ_c’, depending on the CAM diffusivity. When D_CAM is much lower than D_SE (D_CAM/D_SE = 0.001), D_eff rises sharply once a percolation pathway is achieved at the SE fraction above 25%. In contrast, when D_CAM is only one order of magnitude lower than D_SE, lithium ions can diffuse through both materials, as indicated by the non-zero slope before reaching ρ_c’ = 33% in the D_CAM/D_SE = 0.1 case. As a result, there is a non-negligible D_eff before the 33% threshold. For instance, D_eff = 0.25 when the SE fraction is 30%. We hypothesize that above the percolation threshold, the diffusion occurs primarily through the solid electrolyte rather than CAM. We emphasize that the percolation threshold for the SE volume fraction indicates a change in the diffusion behavior, rather than an optimal SE volume fraction for cathode design.

It should be noted that numerical artifacts based on the Fourier-spectral iterative perturbation solver appear when the contrast of material coefficients is large, as indicated in Supplementary Figure 4, where there are fluctuating data below the percolation threshold in the range of D_eff < 10^-4. This issue has been discussed in an analogous study of effective thermal conductivity^[55]. A detailed exploration of these numerical issues is beyond the scope of the present work.

The framework of our microstructure-aware models and the GEMT is general, allowing it to be applied to specific CAM-SE systems, such as NMC-LLZO^[56], as well as sulfide and halide solid electrolytes^[57]. Furthermore, there is growing interest in the development of composite cathodes with gradient phase compositions^[58], and we claim that the ternary diagrams and the analytical model presented in this work can serve as practical tools to accelerate the design and optimization of such systems. With the analytical models, one can quickly generate ternary diagrams with varied percolation thresholds and slopes, enabling the exploration of different transport scenarios. For instance, Supplementary Figure 5 demonstrates examples where the percolation thresholds ρ_c or ρ_c’ are doubled or halved, illustrating how changes in these parameters influence effective transport behavior. These variations, achieved through a unified analytical framework, provide a versatile tool for interpreting transport phenomena across a wide range of three-phase electrochemically active systems, such as solid-state electrolytes with pores and grain boundaries, porous electrodes and separators, and multiphase catalysts. This adaptability makes the analytical models valuable for guiding material design in diverse fields such as energy storage, chemical processing, and environmental remediation. In addition, the flexibility of our model allows us to study a variety of solid electrolytes and active materials. It can also be used to examine effective electronic conductivity based on microstructure with carbon binders/additives.

Beyond transport properties, mechanical considerations also play a critical role in ensuring the stability of composite cathodes^[59]. Our microstructure-resolved approach has the potential to predict effective elastic properties and distribution of stress hotspots, which are key factors in determining the mechanical behavior for solid-state batteries. A similar approach can be applied to solid-state battery materials to address mechanical challenges in cathode design. Specifically, the mesoscale modeling framework can be extended by assigning appropriate elastic constants to each phase and solving the equilibrium equations under applied load. This enables the extraction of effective moduli using homogenization techniques, as well as the identification of local stress concentrations at CAM-SE interfaces and within CAM particles. Such a framework enables high-throughput evaluation of how microstructural features-such as phase fraction, particle shape, and pore morphology-influence both bulk stiffness and the likelihood of mechanical failure, thereby providing design guidelines for mechanically robust solid-state cathodes. This capability has been demonstrated in analogous systems such as metal hydrides encapsulated in host materials^[60].