Phase instability-coupled fracture behavior in garnet LLZO solid electrolytes: a machine learning-enabled atomistic study

Abstract

Fracture in the garnet-type solid electrolyte Li₇La₃Zr₂O₁₂ (LLZO) poses a critical threat to both the performance and safety of solid-state batteries. To unravel the coupled chemomechanical processes that govern fracture in LLZO under external loading, we carry out large-scale molecular-dynamics simulations with a validated machine-learning force field. Our results demonstrate that triaxial stresses at crack flanks trigger a localized cubic-to-tetragonal phase transformation, which is accompanied by Li-ion rearrangement. The emergent tetragonal domains feature lattice contraction normal to the fracture plane, imposing coherent misfit strains that provide an additional driving force for further crack propagation. Crucially, introducing Li deficiencies stabilizes the cubic phase, postponing the phase transition and thereby delaying fracture initiation. These findings highlight the role of intrinsic phase instability in dictating LLZO's fracture resistance and its critical connection to local Li concentration. This chemomechanical coupling points toward targeted strategies to enhance the mechanical robustness of garnet electrolytes, including tuning Li content, ensuring dopant homogeneity, and refining processing protocols.

---

Broader context

The transition to electric vehicles and grid-scale storage relies on batteries that are both energy-dense and intrinsically safe. Solid-state batteries, which replace flammable liquid electrolytes with solid materials, are a promising route but are held back by mechanical failure of the solid electrolyte. A key candidate material, the garnet Li₇La₃Zr₂O₁₂ (LLZO), can crack in service, threatening both lifetime and safety, yet the atomistic origins of this fracture remain unclear. In this work, we use advanced computer simulations to show how local stress and subtle changes in lithium content can trigger a hidden phase change in LLZO that, in turn, drives crack growth. We reveal that making LLZO slightly lithium-deficient stabilizes its structure and delays fracture. These insights provide practical guidelines for industry and researchers on how to tune composition, doping, and processing to improve garnet electrolytes. More broadly, our approach can be extended to study and mitigate failure in other solid-state battery materials.

Introduction

Solid-state batteries (SSBs) are at the forefront of next-generation electrochemical energy storage technologies, promising high energy density and enhanced safety across a wide range of operating conditions.^1–3 The essential component of SSBs is the solid electrolyte (SE), which must combine high thermal and electrochemical stability with fast Li-ion transport, while acting as a physical barrier against Li dendrite growth. Among the many SE chemistries under study,^4–6 the garnet Li₇La₃Zr₂O₁₂ (LLZO) is particularly attractive owing to its high ionic conductivity (in the cubic polymorph) and compatibility with Li metal anodes.^7–9 Yet practical challenges remain, mostly attributable to mechanical degradation in working cells, such as crack formation and propagation under accumulated cycling stresses, coupled with Li dendrite penetration that eventually leads to short circuiting and cell failure.^10–13 Multiple intertwined factors collectively contribute to the mechanical failure of LLZO. For example, the volumetric swing of the cathode due to repeated Li insertion and extraction can exert nonnegligible mechanical stresses on adjacent SEs.¹⁴ Meanwhile, Li diffusion and accompanying redox reactions at LLZO/electrode interfaces can disrupt local bonding, promoting subcritical crack initiation and growth that contributes to cyclic fatigue.^15,16 Applying external stack pressure can help maintain interfacial contact and suppress void formation^17–20 but may elevate mechanical loading on the SE, hastening fracture. Local chemical inhomogeneity, such as excess Li accumulation at grain boundaries (GBs) or pre-existing lattice imperfections, can further impel Li dendrite growth and inter-/trans-granular crack propagation.^21–23 These complications stress the necessity for a fundamental understanding of fracture mechanisms in LLZO, both to forestall short-circuiting and thermal runaway and to design more resilient SEs for durable, high-performance SSBs.^9,10,24,25

To date, most efforts to characterize and predict LLZO's fracture toughness have focused on assessing its inherent brittleness, through either experimental indentation tests^26–31 or Griffith's theory-based modeling.^32–34 To explicitly address the potential impact of GBs, Yu and Siegel constructed atomistic models that suggested significant softening of the elastic response at selected GBs.³⁵ On the experimental front, Fincher et al. combined operando microscopy with dynamic mechanical loading to quantify the residual stresses required to prevent short-circuiting of LLZO.¹³ Athanasiou et al. employed photoelasticity to map out the stress fields around Li dendrites in Ta-doped LLZO during electrochemical cycling, concluding that the stress distributions near Li dendrite tips match those expected to cause brittle fracture of SEs.³⁶ These studies have been very useful in providing plausible microscopic explanations for mechanical degradation modes in LLZO during operation. However, they did not resolve the fundamental atomistic driving forces for fracture initiation, nor did they investigate coupled chemomechanical effects that may critically contribute to crack propagation.

Additional complications arise from the intrinsic metastable nature of cubic LLZO, which is the predominant phase needed for SE applications.³⁷ The conventional cubic phase of LLZO adopts a garnet-type lattice structure (space group Ia3̄d) with random and partial Li⁺ occupation at the tetrahedral (24d) and distorted octahedral (96h) sites that allows for fast Li⁺ diffusion. Upon elemental doping (e.g., Ga and Al) or externally applied stress, LLZO can also crystallize into the acentric cubic phase with the I4̄3d space group.^38,39 In both cubic polymorphs, oxygen anions form a rigid sublattice of corner-sharing ZrO₆ octahedra and LaO₈ dodecahedra, with Li ions threading through large interstitial channels. The high symmetry of the cubic lattice effectively yields nearly isotropic Li⁺ conduction pathways and a delocalized Li⁺ distribution, resulting in high conductivities up to 10⁻⁴–10⁻³ S cm⁻¹ at room temperature.⁴⁰ In contrast, the tetragonal phase (space group I4₁/acd), which is thermodynamically more stable at room temperature, evidences an ordered Li ion arrangement with lowered symmetry: Li⁺ fully populates the 8a, 16f, and 32g sites, leading to a slightly tilted ZrO₆ framework and a broken Li⁺ percolation pathway with substantially reduced ionic conductivity (by two orders of magnitude).

In practice, engineering strategies are often needed to steer the formation of the superionic yet thermodynamically metastable cubic phase over the tetragonal variant. Examples include elemental doping and application of external stress during synthesis.⁴¹ The relative phase stability of LLZO can also be sensitive to Li stoichiometry. For example, local Li depletion or accumulation experienced during repeated cycling or resulting from interfacial side reactions with electrodes can shift local chemical potentials sufficiently to nucleate tetragonal domains.³⁷ Moreover, the off-stoichiometry of Li caused by high-energy ball milling, sintering, or pellet pressing can further promote Li/vacancy ordering,⁴² favoring tetragonal phase formation.

The structural richness of LLZO, with its phase sensitivity to local atomic arrangements of Li and accompanying stress fields, implies that the system could exhibit unique mechanical properties and failure modes during operation. Unfortunately, most prior works have centered on the impact of polymorphic phases on Li-ion transport,^43–46 with only a handful of studies attempting to address the potential coupling with mechanical response. Among these very few reports, Monismith et al. proposed that uniaxial loading-induced cubic-to-tetragonal (c-to-t) phase transformation of LLZO could deteriorate Li-ion transport kinetics^47,48 but act as a toughening mechanism that improves the overall fracture resistance of polycrystalline LLZO.³⁴ In contrast, Hong et al. concluded that mechanical polishing in the presence of excess Li can produce a surface layer of mixed cubic and tetragonal phases that is prone to Li penetration and early-stage short-circuiting.³⁷ In light of this discussion, a clear atomistic elucidation of how the relative stability of LLZO's polymorphs governs fracture resistance and the coupled chemomechanical response is essential. Such an understanding can also rationalize design strategies for enhancing the mechanical robustness of LLZO in practice.

To fill this gap, we carried out large-scale molecular-dynamics simulations using a rigorously benchmarked machine-learning force field (MLMD), aiming to uncover the atomic origins of LLZO's mechanical failure and the role of spontaneous phase evolution during crack initiation and propagation. We considered three LLZO variants: stoichiometric cubic (c-LLZO, Li₇La₃Zr₂O₁₂), Li-deficient cubic (Li_6.3La₃Zr₂O₁₂), and tetragonal (t-LLZO), to address how Li stoichiometry and phase structure influence the fracture behavior of LLZO. By closely analyzing the dynamic evolution of global stress, local lattice parameters and stress–strain fields, as well as local Li ion rearrangement during fracture, we elucidate how stress-induced c-to-t phase transformations at crack flanks contribute to the driving force for crack propagation, thereby compromising the overall fracture resistance of LLZO. Furthermore, we demonstrate the role of Li concentration in dictating LLZO's fracture behavior and propose practical mitigation strategies to stabilize the cubic phase of LLZO for enhanced fracture resistance and cycling longevity.

Results

To allow systematic comparison and elucidation of the impact of phase and composition variations in LLZO, c-LLZO, c-Li_6.3La₃Zr₂O₁₂, and t-LLZO specimens were prepared with the same geometry and initial crack flaw configuration and were subjected to the same externally applied strain field. Prior to the direct fracture simulations on the three LLZO variants, benchmark tests on c-LLZO samples pre-cracked along three orientations (010)[001], (1̄10)[110], and (01̄1)[111] were first conducted to determine the preferential cleavage plane (Fig. S1; SI Movie SI). These dynamical simulations suggest that intragranular fracture in LLZO proceeds via brittle cleavage, with a strong preference for crack opening along the {100} plane. Guided by this result, the three LLZO variants were then initialized with a centered (010)[001] crack of identical dimensions and composition (Fig. 1a). A static volumetric expansion of 6.4%, which is slightly below the yield point of c-LLZO and t-LLZO, determined from hydrostatic expansion tests on their pristine samples (Fig. 1b), was imposed to drive steady-state crack growth while eliminating rate effects. Note that the imposed far-field loading is not intended as a direct proxy for a macroscopic “operating stress” in a full cell. Rather, the mechanism studied here is governed by the local stress state near sharp stress concentrators at nanometer length scales, where strong stress amplification can produce multi-GPa stresses within the crack-tip process zone even under comparatively modest macroscopic loading. Periodic boundary conditions were applied in all directions, and more details regarding the simulation setup can be found in the Methodology section.

Fig. 1 (a) Geometries of the three LLZO variants, i.e., the cubic (c-LLZO and c-Li_6.3La₃Zr₂O₁₂) or tetragonal LLZO (t-LLZO). All initiated with a (001)[001] center crack of the same dimension. (b) Hydrostatic tension–volumetric strain curves of the pristine c-LLZO and t-LLZO samples under isotropic expansion. The dashed lines indicate the volumetric strain (6.4%) applied to the cracked LLZO samples in panels (c1–c3). (c1–c3) Stress-release curves of c-LLZO, t-LLZO, and c-Li_6.3La₃Zr₂O₁₂ in the [100], [001], and [010] directions during fracture simulations.

Fig. 1(c1–c3) show the global stress-release responses at 300 K for the three LLZO variants along [100], [010], and [001]. Interestingly, the pristine samples behave nearly identically up to 6.4% volumetric strain, but introducing a crack flaw produces marked divergence. The tetragonal phase (t-LLZO; blue curve) shows no stress drop over the entire simulation. By contrast, cubic LLZO (c-LLZO; yellow curve) exhibits a two-stage response along [100] and [001]: stress plateaus for ∼20 ps (stage I) and then declines (stage II), whereas along [010] it decays monotonically. The c-Li_6.3La₃Zr₂O₁₂ composition (red curve) likewise shows a two-stage response, with a pronounced in-plane stress buildup during stage I followed by release in stage II. Although c-Li_6.3La₃Zr₂O₁₂ appears to absorb less mechanical work under hydrostatic expansion (smaller area under the curve in Fig. 1b), its crack propagation is delayed relative to c-LLZO (Fig. 1c), suggesting improved mechanical stability in the presence of pre-existing cracks.

These differences in mechanical response among the three LLZO variants are also evidenced in atomistic renderings of bond rupture, as shown in Fig. 2 (focusing on the Zr–O bonding network for illustration purposes). In pristine LLZO, each Zr⁴⁺ center is octahedrally coordinated by six O²⁻ ions, forming a rigid 3-dimensional, interconnected backbone structure. During fracture of c-LLZO (Fig. 2a), ZrO₆ octahedra at the crack surface progressively convert into ZrO₅ and then ZrO₄, with the bond-breakage front propagating from the crack sides toward the tips. In contrast, t-LLZO (Fig. 2b) retains the ZrO₆ coordination environment throughout, underscoring its superior lattice resilience to crack propagation. The Li-deficient variant (Fig. 2c) also undergoes a Zr–O bond rupture, but the rate at which ZrO₆ octahedra break down is markedly slower than that in c-LLZO. This slower bond-rupture progression aligns with the delayed stress-release behavior observed in Fig. 1(c1–c3).

Fig. 2 Illustrations of Zr–O bonding states around the crack in (a) c-LLZO, (b) t-LLZO, and (c) c-Li_6.3La₃Zr₂O₁₂.

Since no plane slip, glide, or dislocation activity was observed in the Zr–O (or La–O) sublattices of any variant, we propose that the deteriorated fracture resistance of c-LLZO and c-Li_6.3La₃Zr₂O₁₂ relative to t-LLZO arises from transient disorder in the Li sublattice and its stress-driven dynamics. Within this interpretation, the stage-I stress stagnation in Fig. 1(c1 and c2) reflects the ionic anelasticity of Li in the cubic polymorphs—that is, a time-dependent but largely recoverable strain produced by stress-induced redistribution/hopping of mobile Li among metastable sites. This relaxation is thermally activated and operates on picosecond–nanosecond timescales relevant to our simulations. By contrast, the more ordered Li arrangement in t-LLZO suppresses such anelastic relaxation, so its mechanical response remains predominantly elastic and the local strength of the Zr–O network governs crack-opening resistance. Consistent with the thermal origin of Li mobility and anelasticity, a comparable c-LLZO fracture simulation at 0 K sustains high stress without crack advance (green curves in Fig. 1(c1–c3)), confirming that thermally activated processes degrade fracture resistance in the cubic polymorphs.

Nevertheless, Li mobility or near-crack bond rupture alone cannot explain the anisotropic stress evolution observed in stage I of the cubic polymorphs in Fig. 1(c1–c3). To uncover the additional mechanisms at play, we performed a more detailed analysis of the local lattice structure, atomic strain, and virial stress evolution during fracture. As shown in Fig. 3a, adaptive common neighbor analysis (CNA)^49,50 of the Zr sublattice reveals that a body-centered symmetry (BCC or BCT) persists throughout crack propagation. Taking advantage of this well-defined Zr sublattice, we computed the local lattice “c/a” ratio from the atomic elastic strain tensor, π_c/a = (1 + ε_yy)/(1 + ε_xx) (see SI S-I for calculation details) as a proxy to quantify local lattice distortion. As shown in Fig. 3b, regions where π_c/a > 1 appear at the crack tip and its adjacent wing-shaped regions, whereas π_c/a < 1 emerges along the crack sides due to local biaxial tension. One might initially map any deviation from π_c/a = 1 to a tetragonal domain (π_c/a = 0.96 for t-LLZO), but it remains inconclusive whether these distortions arise solely from elastic deformation of the cubic lattice or are indeed associated with the phase evolution from cubic to tetragonal.

Fig. 3 Evolution of the local atomic structure in c-LLZO during fracture. (a) Zr sublattice symmetry from adaptive common-neighbor analysis; body-centered cubic (BCC) or body-centered tetragonal (BCT) environments are shown in blue, while sites lacking clear symmetry (“frustrated” structures) are rendered in gray. (b) Spatial map of the local lattice c/a ratio, computed from the Zr sublattice's local elastic strain (a∥x, c∥y). (c) Local volumetric strain of the Zr-sublattice; pink markers denote regions of negative (compressive) volumetric change. (d) Distribution of local atomic viral stress. In each color bar, the white line marks the zero or middle value.

To distinguish between elastic strain and genuine phase change, we also computed local volumetric strain from instantaneous lattice parameters. Surprisingly, negative volumetric strain pockets appear along the crack flanks (highlighted by the pink spots in Fig. 3c). Following a pure linear-elastic response theory under tensile loading, one would expect positive volume dilation, and therefore these localized volume contractions indicate a potential phase evolution. Further examination of the local atomic stress distribution (σ_yy) in c-LLZO (Fig. 3d; see Fig. S2 for σ_xx and σ_xy) highlights that tensile stresses concentrate at the crack tips, effectively driving smooth crack propagation along the (001) plane. At the same time, pronounced compressive stress regions flank the fracture plane, co-located with the negative volumetric strain pockets. This spatial overlap confirms a local c-to-t phase transformation of LLZO that is directly accompanied by crack propagation.

Having identified the local c-to-t transformation, we next verified that the stress-induced tetragonal phase corresponds to the well-established I4₁/acd space group and examined how Li-ion stoichiometry influences this phase transformation behavior. As noted earlier, Li ions adopt a more ordered arrangement in t-LLZO than in c-LLZO. To quantify this ordering, we employed a configurational entropy fingerprint⁵¹ that projects integrals of the Li–Li radial distribution functions onto individual Li atoms (see SI S-II for calculation details). The resulting entropy values are always negative, with more negative values indicating greater local order/crystallinity and less negative values indicating more disordered/liquid-like atomic distributions. During the c-to-t phase transformation, a decrease in the local entropy is thus expected for the disorder-to-order Li-ion rearrangement (see Fig. S3 for statistical Li distributions in each phase). Fig. 4(a1 and a2) reveal that as a crack propagates in c-LLZO, Li ions along the crack flanks undergo a pronounced decrease in entropy, reflecting their rearrangement into a more ordered, tetragonal-like lattice structure. A milder entropy reduction appears in the wing-shaped regions ahead of the crack tips, whereas the immediate tip zone exhibits higher entropy, reflecting local structural frustration. In the Li-deficient variant c-Li_6.3La₃Zr₂O₁₂ (Fig. 4(a3 and a4)), the same disorder-to-order transition occurs, but with a clear delay relative to c-LLZO. Quantitative analysis of Li-site entropy at crack tips and sides further confirms that c-Li_6.3La₃Zr₂O₁₂ lags behind c-LLZO in this Li-rearrangement process (Fig. S4), consistent with the observed prolonged stage I stress plateau in Fig. 1(c1–c3).

Fig. 4 Li-ion rearrangement in initially cubic phase LLZO during fracture. (a1–a4) Evolution of local entropy for Li ions in c-LLZO and c-Li_6.3La₃Zr₂O₁₂, respectively; only Li sites are displayed, with cooler colors indicating lower (more negative) entropy and hence greater ordering. (b) Site-specific entropy signatures in pristine t-LLZO: top row shows per-atom entropy for Li at the three Wyckoff positions (8a, 16f, and 32g); middle row displays Li with entropy values below −10 in t-LLZO (aligned atom sites denote 8a-site Li in body-centered symmetry); bottom row shows pristine c-LLZO using the same threshold, with atom sites lacking identifiable symmetry. (c1 and c2) Li ions with local entropy values below −10 in c-LLZO and c-Li_6.3La₃Zr₂O₁₂ during fracture; blue atoms denote 8a-site Li in body-centered symmetry, grey atoms denote Li sites lacking identifiable symmetry and red outlines highlight the advancing crack fronts.

To further resolve Li ordering at specific atomic sites, we deconvoluted the entropy metric on the well-defined high-symmetry Wyckoff positions (8a, 16f, and 32g) within the t-LLZO lattice setting. As shown in Fig. 4b (top row), each site exhibits a unique entropy signature, with Li atoms at the 8a positions showing the lowest entropy values. Setting an entropy threshold excludes the 8a-site Li contribution, thereby highlighting the body-centered signature that is unique to the t-LLZO symmetry (Fig. 4b middle row) and absent for the cubic symmetry (Fig. 4b bottom row). Applying the same threshold to the fracture simulation results reveals a body-centered Li arrangement along the crack sides in c-LLZO as the crack extends (Fig. 4c1). Similar behavior is also evident in c-Li_6.3La₃Zr₂O₁₂, although the onset is delayed (Fig. 4c2; see SI Movie SII for complete dynamics). Notably, both stoichiometric and Li-deficient cubic LLZO preserve isolated regions lacking body-centered order, suggesting that the transformed tetragonal patches retain a certain degree of Li sublattice disorder. Accordingly, even in the Li-deficient sample, Li re-ordering is observed within individual LLZO unit cells. This allows the c-to-t transformation to proceed on tens-of-picoseconds timescales, distinct from typical diffusive phase transformations that involve slow, long-range ion migration. We thus propose that this rapid, local Li rearrangement before crack propagation is primarily responsible for the stage-I stress plateau in Fig. 1(c1–c3). Although not explicitly modeled here, we acknowledge that Li redistribution on longer diffusive timescales in the vicinity of crack-like flaws could further modify local Li ordering and, consequently, the propensity for phase transformation. We also note that finite-size, finite-time atomistic simulations cannot fully capture slow nucleation pathways that may operate at much lower driving forces or across microstructural length scales.

We now proceed to discuss how Li concentration may affect the phase transformation behavior of LLZO and revise the stress distribution profile. Fig. 5a plots the calculated absolute free energy⁵² of both stoichiometric and Li-deficient LLZO at 300 K against the c/a ratio along the Bain path. For stoichiometric LLZO, the tetragonal polymorph (c/a = 0.96) lies at the global energy minimum, whereas the metastable (dashed parabolas) cubic configurations reside at higher-energy states. The energy barrier for c-to-t transformation is relatively small (∼2 kJ mol⁻¹), whereas the backward t-to-c transformation requires significantly higher energy (∼40 kJ mol⁻¹). Such an asymmetric energy landscape implies that modest axial stresses can readily trigger c-to-t transformations during fracture. By comparison, Li-deficient c-Li_6.3La₃Zr₂O₁₂ exhibits the cubic polymorph as its lowest-energy phase, confirming that Li vacancies can indeed stabilize the cubic lattice and thereby enhance resistance to the stress-induced c-to-t phase transition.

Fig. 5 Influence of Li stoichiometry on LLZO phase stability and fracture response. (a) Free-energy profiles along the Bain path at 300 K for stoichiometric and Li-deficient LLZO, each plotted relative to its own minimum. (b) Uniaxial [001] stress–strain curves under tension and compression for pristine c-LLZO, t-LLZO, and c-Li_6.3La₃Zr₂O₁₂. (c) Schematic of crack advancement in c-LLZO if the cubic phase were to remain stable (left), in comparison with the observed scenario when a local c-to-t phase transformation occurs at the crack flanks (second from left); arrows indicate the additional driving force opening the crack from coherent misfit strains.

Motivated by these free-energy analyses, we carried out additional direct rate-dependent uniaxial tensile loading simulations on the three LLZO variants to confirm the stress-driven phase evolution behavior of LLZO. In contrast to their very similar response under hydrostatic tension (Fig. 1b), the stress–strain curves obtained for pristine c-LLZO, t-LLZO, and c-Li_6.3La₃Zr₂O₁₂ under uniaxial loading diverge markedly (Fig. 5b). In c-LLZO under compression and in t-LLZO under tension, pronounced stress plateaus (indicated by arrows in Fig. 5b) imply energy dissipation via a phase-transition mechanism. Moreover, the c-LLZO curve deviates from perfect linear elasticity at tensile strains as low as ∼2%, likely related to early activation of Li diffusion between the interstitial sites and incipient c-to-t lattice distortions inherent to the metastable cubic polymorph. Again, Li deficiency in c-Li_6.3La₃Zr₂O₁₂ suppresses the c-to-t phase transformation under axial loading, confirming that the reduced Li content helps stabilize the cubic phase against applied stress. These orientation-dependent mechanical responses explain why the cubic polymorph exhibits anisotropic stress evolution in stage I (Fig. 1(c1–c3)). Under ideal, flaw-free conditions, LLZO experiences uniform hydrostatic stresses, which are not likely to drive a c-to-t transformation. However, in an operating battery with a cathode composite, LLZO grains rarely experience truly hydrostatic stress. In this case, cracks introduce free surfaces that shift the local stress state from (near-)hydrostatic to triaxial. Such triaxial stresses are capable of triggering the c-to-t transformation, accounting for the anisotropic stress response seen in stage I.

Discussion

Based on our simulation results, we identify three unique features that affect LLZO's fracture behavior. First, the highly mobile Li ion sublattice in cubic LLZO introduces anelastic, thermally activated processes that cause deviation from a linear-elastic fracture mechanism. According to Griffith's criterion,⁵³ the fracture toughness is given by , where γ is the surface energy, E is the directional Young's modulus, and ν is Poisson's ratio. Using parameters extracted from our MLMD simulations (Table S1), this estimation yields K_IC = 0.525 MPa m^−0.5 for c-LLZO, 0.522 MPa m^−0.5 for c-Li_6.3La₃Zr₂O₁₂, and 0.532 MPa m^−0.5 for t-LLZO. In contrast, our dynamic MLMD simulations rank fracture resistance as t-LLZO ≫ c-Li_6.3La₃Zr₂O₁₂ > c-LLZO, underscoring the insufficiency of the linear-elastic assumption for LLZO.

Second, tension–compression asymmetry arising from stress-induced phase transformations significantly alters local stress/strain fields around the crack flaws. Fig. 5c schematizes how the local c-to-t phase transformation degrades LLZO's fracture resistance. In the absence of phase evolution, the fracture response of c-LLZO is governed solely by its intrinsic elastic (brittle-cleavage) properties, and the crack advances following the applied stress field and the material's inherent cleavage toughness. However, (bi)axial tensile stress at the crack flanks is found to be sufficient to trigger a localized c-to-t phase transformation, with the newly formed tetragonal phase possessing a reduced lattice constant normal to the fracture plane. This results in coherent misfit strain that induces additional tensile stress at the crack tip, effectively adding to the driving force for crack propagation in a positive feedback cycle. This process stands in stark contrast to transformation-toughened mechanisms, where stress-induced phase changes at the crack tip absorb strain energy and impede crack growth.

Third, the propensity for c-to-t phase evolution under mechanical loading is highly sensitive to local Li occupancy because stoichiometry directly controls the stability of the cubic polymorph. Li deficiency stabilizes the cubic phase and thus tends to suppress stress-assisted phase transformations and improve fracture resistance by reducing transformation-induced misfit strain. In practice, however, Li stoichiometry is unlikely to remain spatially uniform near defects. Hydrostatic tension in the crack-tip stress field can bias Li/vacancy redistribution, so a crack-localized deviation in Li content—toward either enrichment or depletion, depending on boundary conditions and defect chemistry—will shift the local distance to the cubic–tetragonal boundary and thereby modify the transformation onset.

These findings have broader implications for practical applications of LLZO. For instance, elemental doping has been recognized as an effective strategy for stabilizing the cubic phase of LLZO.^44,54–56 However, processing steps such as sintering and volatilization compensation can generate spatially inhomogeneous dopant distributions in the bulk electrolyte,^37,57 leaving locally under-doped domains that sit closer to the cubic–tetragonal boundary and are therefore more susceptible to stress-assisted c-to-t transformation in the triaxial tensile fields surrounding pores, microcracks, and sharp flaws. At interfaces, additional departures from the bulk Li occupancy can arise from space-charge effects, dopant segregation, and interfacial reaction products, producing Li-nonstoichiometric zones whose magnitude and even sign depend on the specific interfacial chemistry and polarization state. Such local phase/stoichiometry variations introduce eigenstrains (misfit strains) and modulus contrasts that can nucleate microcracks or intensify crack-tip driving forces, thereby degrading the effective fracture resistance. Likewise, because higher Li content favors the tetragonal phase,^55,58 current-density hotspots during cycling can induce local Li enrichment, promoting a c-to-t transformation that amplifies stress accumulation surrounding defects.

To counteract phase-instability-coupled mechanical degradation of LLZO, it is useful to first translate the transformation-assisted fracture mechanism identified here into processing-tunable microstructural descriptors. Three descriptors are therefore particularly actionable: (i) porosity and microcrack density, which determine the density/acuity of stress concentrators and thus the likelihood of reaching the triaxial stress threshold; (ii) grain size and texture, which control orientation-dependent stress partitioning and transformation-variant mismatch across GBs and triple junctions; and (iii) dopant/vacancy homogeneity, chemical corrosion,^59,60 and GB segregation, which set local cubic stability by creating under-doped or vacancy-inhomogeneous domains. Therefore, the mitigation strategies should consider minimizing flaw populations and maximizing chemical/defect uniformity. Particularly, stress engineering,⁶¹ crystallographic texturing or coarse-grained microstructures approaching the single-crystal limit,^40,62 and applying amorphous coatings,^63–65 can help suppress these effects by narrowing orientation scatter and reducing the density of high-acuity stress concentrators.

Conclusions

In this study, we carry out a detailed atomic-scale investigation of intragranular fracture in LLZO under applied mechanical loads. Although c-LLZO appears to fail by conventional brittle cleavage, our simulations expose underlying microstructural evolutions that govern its true fracture behavior. We show that the highly triaxial stress state at crack flanks can trigger a local c-to-t transformation that is accompanied by significant Li-ion rearrangement. This rearrangement shrinks the lattice normal to the crack plane, imposing coherent misfit strains that add to the driving force for crack propagation in a positive feedback loop. By varying Li content, we demonstrate that reducing Li concentration delays the onset of c-to-t transformation and correspondingly mitigates the material's susceptibility to fracture. In contrast to conventional mechanical studies that focus on bulk, uniform deformation, our work unveils a novel, phase-transformation-driven mechanism of internal stress generation that critically undermines LLZO's mechanical integrity. These insights suggest practical design guidelines for LLZO-based SEs that are grounded in chemomechanical coupling, such as prioritizing homogeneous dopant and Li distributions to stabilize the cubic phase and mitigating Li-enrichment hotspots during cycling. Together, such strategies promise to enhance both the fracture toughness and ionic-transport performance of LLZO in SSB applications.

Methods

Development and validation of ML interatomic potential for LLZO

An ML interatomic potential was developed following the same framework reported in our previous work⁶⁶ using the same parameters for the neural network architecture (2 hidden layers and 20 nodes in each layer), optimization algorithm, and Behler–Parrinello radial and angular symmetric functions. To enable simulations for mechanical response and fracture propagation, the training dataset was expanded by adding surface data to the existing dataset (5466 structures). The surface data were generated by sampling atomic structures and their total energies and atomic forces every 50 fs or longer from ab initio MD (AIMD) simulations, using the (001) plane with Li termination and the (111) plane with Zr termination (see SI S-III for details of density functional theory (DFT) and AIMD simulation methods). The number of new data points is 1540 for the (001) plane at temperatures between 300 and 3000 K and 975 for the (111) plane at temperatures between 800 and 1600 K, and thus in total 7981 structures exist in the new dataset. The dataset explicitly includes surface disordering and amorphization with a significant decrease in density, which can occur during the process of crack initiation and propagation.

We validated our ML model using energy and force errors for cubic and amorphous structures as well as six slab models of LLZO (the amorphous and six surface structures were adopted from our previous report).⁴⁶ It was done by performing MD simulations with the developed ML potential (in short, MLMD) for 100 ps and comparing the total energies and atomic forces between DFT calculations and ML predictions using atomic structures sampled from MD trajectories every 1 ps. Fig. S5 shows that the root mean squared errors (RMSEs) of energies and forces at high temperatures are less than 1.6 meV per atom and 0.11 eV Å⁻¹ for cubic LLZO at 2000 K and 1800 K with NVT and NPT ensembles, respectively; and 4.6 meV per atom and 0.17 eV Å⁻¹ for amorphous LLZO at 3000 K and 2000 K with NVT and NPT ensembles, respectively. The RMSEs for the six slab models are also small, being less than 8.4 meV per atom and 0.19 eV Å⁻¹ (Fig. S5). Two important points are noted: (1) the ML potential can perform reliable MD simulations for various surface types that it did not learn, e.g., surfaces with the (110) plane with three different terminations and the (111) plane with Li termination as shown in Fig. S6, implying the transferability to simulate various surfaces; and (2) the ML potential can predict a wide range of densities of amorphous LLZO as shown with MD simulations for the slab models at 2000 K in Fig. S6, where LLZO is completely disordered/melted and fills the entire vacuum area with a significant decrease in density (up to a 73.7% increase in volume in our tests), which can happen at the initial stage of crack opening and fracture surfaces with strains. Thus, the ML potential can predict accurate energies and forces of cubic and amorphous LLZO as well as LLZO surfaces with various densities at high temperatures.

We further validated our ML potential by predicting structural, vibrational, dynamical, and elastic properties of LLZO; see our previous report for property calculation methods.⁶⁶ Fig. S7 presents the radial distribution functions (RDFs) of all element pairs for cubic LLZO at 1273 K with +3% hydrostatic strain and amorphous LLZO at 3000 K with +10% strain, showing that MLMD simulations predict nearly the same RDFs compared to AIMD, whereas CMD simulations lead to different RDFs from AIMD. In addition, the ML potential can predict an accurate phase transformation temperature between tetragonal and cubic phases near 900 K (Fig. S8), which is close to the experimental measurement of 918 K,⁶⁷ showing that our ML potential can predict structural characteristics accurately. The predicted vibrational density of states of elements and Li diffusivities in cubic and amorphous LLZO from MLMD simulations are also close to those calculated from AIMD trajectories, whereas CMD simulations result in different characteristics of vibrational and dynamic properties from AIMD (Fig. S9 and S10). Moreover, the predicted bulk, shear, and Young's moduli as well as Poisson's ratio by ML potential are 102.6 GPa, 60.0 GPa, 150.6 GPa, and 0.26, respectively, which are close to those calculated by DFT (114.7 GPa, 65.2 GPa, 164.5 GPa, and 0.26)⁶⁶ as well as experimental measurements (102.8 ± 0.3 GPa, 59.6 ± 0.1 GPa, 149.8 ± 0.4 GPa, and 0.257 ± 0.002 by resonant ultrasound spectroscopy for Li_6.24La₃Zr₂Al_0.24O_11.98 at 298 K).²⁸ The room-temperature power spectra of LLZO from MLMD (Fig. S12) also agree well with the experimental Raman frequencies.⁶⁸ Thus, we conclude that the ML potential can predict the structural, vibrational, dynamical, and elastic properties of cubic and amorphous LLZO accurately.

To evaluate the performance of our developed ML potential for fracture simulations at larger length scales, we compared the mechanical responses of cubic LLZO with a surface edge crack under constant tensile strain, simulated using MLMD and MD with a classical interatomic potential,⁶⁹ respectively. Fig. S11 and SI Movie SIII show that our MLMD simulations yield a reasonably sharp and smooth crack propagation from the surface into the LLZO grain interior along the [100] direction and the (010) surface plane. In contrast, the same system simulated with the classical potential undergoes abnormal surface structural reconfigurations, and within the simulation time frame, the applied strain could not drive crack propagation due to crack tip blunting, which is not realistic. The overall validation and performance tests indicate the superior capability of our MLMD in simulating the fracture behavior of LLZO with high fidelity.

Machine-learning molecular dynamics simulations

Preparation of pristine LLZO samples. Initial unit-cell geometries for both cubic and tetragonal LLZO were taken from our prior DFT calculations.⁶⁶ Each unit cell was replicated to produce supercells of the target dimensions for large-scale fracture MLMD simulations using the LAMMPS package⁷⁰ and our validated ML interatomic potential. All models were first relaxed at 0 K via conjugate-gradient energy minimization, then equilibrated at 300 K and zero isotropic pressure in the NPT ensemble (Nose–Hoover thermostat and barostat) for 40 ps. Periodic boundary conditions (PBCs) were applied in all directions, and a 0.5 fs timestep was used throughout. At 300 K, the relaxed lattice constants for the cubic phase were a = b = c = 13.112 Å, and for the tetragonal phase, a = b = 13.287 Å and c = 12.751 Å.

Crack introduction and fracture simulations. To create identical pre-cracked specimens, we removed all atoms within a central rectangular prism (2 nm long × 0.2 nm wide) on the x–y plane, carefully preserving overall stoichiometry and removing the same set of atomic species in c-LLZO and t-LLZO. These samples were then thermalized at 300 K and zero isotropic pressure under the NPT ensemble for 100 ps. The Li-deficient sample, c-Li_6.3La₃Zr₂O₁₂, was generated by randomly removing 10 at% of Li from the relaxed c-LLZO sample. After equilibration, each variant introduced with the crack was uniformly strained to 6.4% volumetric expansion, equivalent to a 2.1% tensile strain along each Cartesian axis. Subsequent fracture simulations were performed in the NVT ensemble at 300 K, allowing steady-state crack propagation.

Uniaxial loading simulations. Uniaxial loading tests were performed on pristine c-LLZO, t-LLZO, and c-Li_6.3La₃Zr₂O₁₂ samples along the [001] direction at a constant engineering strain rate of 4 × 10⁸ s⁻¹. During each run, the sample was held at 300 K in the NPT ensemble, with the barostat maintaining zero pressure in the lateral [100] and [010] directions while imposing continuous strain along [001].

Surface energy calculations. Surface energies for c-LLZO, t-LLZO, and c-Li_6.3La₃Zr₂O₁₂ were computed as γ = (E_surface − E_bulk)/2A, where E_surface is the total energy of an LLZO slab containing two surfaces, E_bulk is the energy of the fully periodic bulk cell, and A is the area of one surface. Slab models were constructed with the ASE package⁷¹ and relaxed at 300 K in the NVE ensemble with a Langevin thermostat for 0.2 ns. Energies were sampled every 50 fs over the final 5 ps of relaxation to yield converged surface-energy values.

Free energy calculations. Absolute free energies for both stoichiometric and Li-deficient LLZO were obtained via nonequilibrium thermodynamic integration following the Frenkel–Ladd approach as proposed by Freitas et al.⁷² Along the Bain path, we systematically varied the c-axis lattice parameter while coupling a and b to relax at each step to minimize the total energy. To improve statistics, each LLZO unit cell was replicated into a 4 × 4 × 4 supercell, resulting in 12 288 atoms for the Li7.0 composition and 11 939 atoms for the Li6.3 composition. The switching time t_switch and equilibration time t_equil were both set to 5 ps. It is worth noting that Li ions in cubic LLZO diffuse at 300 K; therefore, the Einstein-crystal reference lacks a single characteristic frequency. To address this, we computed the mean-squared displacement, (Δr)², over the final 0.5 ps of each equilibration window (t_total = t_equil + 0.5 × t_switch). These (Δr)² values were used to assign an effective spring constant for the Einstein crystal reference, ensuring a consistent integration protocol across all c/a ratios and both polymorphs.

Conflicts of interest

There are no conflicts to declare.

Data availability

The data supporting this article have been included as part of the supplementary information (SI). Supplementary information is available. See DOI: https://doi.org/10.1039/d6eb00053c.

Acknowledgements

This work was sponsored by the Department of Energy Transportation Technologies Office, and was performed under the auspices of the U.S. Department of Energy by Lawrence Livermore National Laboratory under contract DE-AC52-07NA27344. Computing support was provided by the Department of Energy, located at the National Laboratory of the Rockies and the Computing Grand Challenge program from Lawrence Livermore National Laboratory.

References

1. T. Thompson, S. H. Yu, L. Williams, R. D. Schmidt, R. Garcia-Mendez, J. Wolfenstine, J. L. Allen, E. Kioupakis, D. J. Siegel and J. Sakamoto, ACS Energy Lett., 2017, 2, 462–468.
2. S. Xin, Y. You, S. Wang, H. C. Gao, Y. X. Yin and Y. G. Guo, ACS Energy Lett., 2017, 2, 1385–1394.
3. Y. Y. Ren, T. Danner, A. Moy, M. Finsterbusch, T. Hamann, J. Dippell, T. Fuchs, M. Müller, R. Hoft, A. Weber, L. A. Curtiss, P. Zapol, M. Klenk, A. T. Ngo, P. Barai, B. C. Wood, R. P. Shi, L. W. F. Wan, T. W. Heo, M. Engels, J. Nanda, F. H. Richter, A. Latz, V. Srinivasan, J. Janek, J. Sakamoto, E. D. Wachsman and D. Fattakhova-Rohlfing, Adv. Energy Mater., 2023, 13, 2201939.
4. S. Ramakumar, C. Deviannapoorani, L. Dhivya, L. S. Shankar and R. Murugan, Prog. Mater. Sci., 2017, 88, 325–411.
5. K. Kerman, A. Luntz, V. Viswanathan, Y. M. Chiang and Z. B. Chen, J. Electrochem. Soc., 2017, 164, A1731–A1744.
6. Z. Z. Zhang, Y. J. Shao, B. Lotsch, Y. S. Hu, H. Li, J. Janek, L. F. Nazar, C. W. Nan, J. Maier, M. Armand and L. Q. Chen, Energy Environ. Sci., 2018, 11, 1945–1976.
7. Z. D. Gordon, T. Yang, G. B. Gomes Morgado and C. K. Chan, ACS Sustainable Chem. Eng., 2016, 4, 6391–6398.
8. E. A. Il'ina, A. A. Raskovalov and O. G. Reznitskikh, J. Chem. Thermodyn., 2019, 128, 68–73.
9. Y. Wang, Z. Chen, K. Jiang, Z. X. Shen, S. Passerini and M. H. Chen, Small, 2024, 20, 2402035.
10. L. Porz, T. Swamy, B. W. Sheldon, D. Rettenwander, T. Frömling, H. L. Thaman, S. Berendts, R. Uecker, W. C. Carter and Y. M. Chiang, Adv. Energy Mater., 2017, 7, 1701003.
11. H. K. Tian, B. Xu and Y. Qi, J. Power Sources, 2018, 392, 79–86.
12. X. M. Liu, R. Garcia-Mendez, A. R. Lupini, Y. Q. Cheng, Z. D. Hood, F. D. Han, A. Sharafi, J. C. Idrobo, N. J. Dudney, C. S. Wang, C. Ma, J. Sakamoto and M. F. Chi, Nat. Mater., 2021, 20, 1485–1490.
13. C. D. Fincher, C. E. Athanasiou, C. Gilgenbach, M. Wang, B. W. Sheldon, W. C. Carter and Y. M. Chiang, Joule, 2022, 6, 2794–2809.
14. R. Mücke, M. Finsterbusch, P. Kaghazchi, D. Fattakhova-Rohlfing and O. Guillon, J. Power Sources, 2021, 489, 229430.
15. M. Klinsmann, F. E. Hildebrand, M. Ganser and R. M. McMeeking, J. Power Sources, 2019, 442, 227226.
16. L. Barroso-Luque, Q. S. Tu and G. Ceder, J. Electrochem. Soc., 2020, 167, 020534.
17. M. J. Wang, R. Choudhury and J. Sakamoto, Joule, 2019, 3, 2165–2178.
18. C. H. Yuan, X. Gao, Y. K. Jia, W. Zhang, Q. L. Wu and J. Xu, Nano Energy, 2021, 86, 106057.
19. Q. S. Tu, L. Barroso-Luque, T. Shi and G. Ceder, Cell Rep. Phys. Sci., 2020, 1, 100106.
20. J. M. Doux, H. Nguyen, D. H. S. Tan, A. Banerjee, X. F. Wang, E. A. Wu, C. Jo, H. D. Yang and Y. S. Meng, Adv. Energy Mater., 2020, 10, 1903253.
21. H. K. Tian, Z. Liu, Y. Z. Ji, L. Q. Chen and Y. Qi, Chem. Mater., 2019, 31, 7351–7359.
22. C. Zhu, T. Fuchs, S. A. L. Weber, F. H. Richter, G. Glasser, F. Weber, H. J. Butt, J. Janek and R. Berger, Nat. Commun., 2023, 14, 1300.
23. S. Heo, D. Lee, K. Kim, Y. Kim, D. J. Yun, S. Park, J. Lee, S. Kim, J. S. Kim and S. Park, J. Power Sources, 2021, 510, 230389.
24. B. C. Wood, J. B. Varley, K. E. Kweon, P. Shea, A. T. Hall, A. Grieder, V. P. Aguirre, D. Rigling, E. L. Ventura, C. Stancill and N. Adelstein, Philos. Trans. R. Soc., A, 2021, 379, 20190467.
25. S. Kalnaus, N. J. Dudney, A. S. Westover, E. Herbert and S. Hackney, Science, 2023, 381, eabg5998.
26. Z. Z. Fu, L. Zhang, J. E. Gritton, G. Godbey, T. Hamann, Y. H. Gong, D. McOwen and E. Wachsman, ACS Appl. Mater. Interfaces, 2020, 12, 24693–24700.
27. Y. Kim, H. Jo, J. L. Allen, H. Choe, J. Wolfenstine and J. Sakamoto, J. Am. Ceram. Soc., 2016, 99, 1367–1374.
28. J. E. Ni, E. D. Case, J. S. Sakamoto, E. Rangasamy and J. B. Wolfenstine, J. Mater. Sci., 2012, 47, 7978–7985.
29. J. Wolfenstine, H. Jo, Y. H. Cho, I. N. David, P. Askeland, E. D. Case, H. Kim, H. Choe and J. Sakamoto, Mater. Lett., 2013, 96, 117–120.
30. F. Aguesse, W. Manalastas, L. Buannic, J. M. L. del Amo, G. Singh, A. Llordés and J. Kilner, ACS Appl. Mater. Interfaces, 2017, 9, 3808–3816.
31. J. Zhao, Y. F. Tang, Q. S. Dai, C. C. Du, Y. Zhang, D. C. Xue, T. W. Chen, J. Z. Chen, B. Wang, J. M. Yao, N. Zhao, Y. S. Li, S. M. Xia, X. X. Guo, S. J. Harris, L. Q. Zhang, S. L. Zhang, T. Zhu and J. Y. Huang, Energy Environ. Mater., 2022, 5, 524–532.
32. S. Monismith, C. D. Fincher, Y. M. Chiang, J. Qu and R. Dingreville, Adv. Energy Mater., 2024, 14, 2303567.
33. J. Wolfenstine, J. L. Allen, J. Sakamoto, D. J. Siegel and H. Choe, Ionics, 2018, 24, 1271–1276.
34. S. Monismith, J. Qu and R. Dingreville, J. Mech. Phys. Solids, 2022, 160, 104791.
35. S. Yu and D. J. Siegel, ACS Appl. Mater. Interfaces, 2018, 10, 38151–38158.
36. C. E. Athanasiou, C. D. Fincher, C. Gilgenbach, H. J. Gao, W. C. Carter, Y. M. Chiang and B. W. Sheldon, Matter, 2024, 7, 95–106.
37. S. K. Hong, K. H. Shin, S. Kim, S. H. Song, K. S. Kim, D. Lee, S. H. Yu, S. K. Jung and H. Kim, Adv. Energy Mater., 2025, 15, 2402666.
38. X. Xiang, Z. C. Fang, F. Chen, H. H. Wang, W. Y. Yang, C. H. A. Wei, J. B. Yang, X. B. Ma, D. F. Chen, K. Sun, Y. H. Zhang and Q. Shen, Ceram. Int., 2022, 48, 9371–9377.
39. M. B. Dixit, B. S. Vishugopi, W. Zaman, P. Kenesei, J. S. Park, J. Almer, P. P. Mukherjee and K. B. Hatzell, Nat. Mater., 2022, 21, 1298–1305.
40. M. Philipp, B. Gadermaier, P. Posch, I. Hanzu, S. Ganschow, M. Meven, D. Rettenwander, G. J. Redhammer and H. M. R. Wilkening, Adv. Mater. Interfaces, 2020, 7, 2000450.
41. M. M. Raju, F. Altayran, M. Johnson, D. Wang and Q. Zhang, Electrochem, 2021, 2, 390–414.
42. Z. H. Yan and Y. Z. Zhu, Chem. Mater., 2024, 36, 11551–11557.
43. K. Meier, T. Laino and A. Curioni, J. Phys. Chem. C, 2014, 118, 6668–6679.
44. X. Cheng, J. Huang, W. Qiang and B. Huang, Ceram. Int., 2020, 46, 3731–3737.
45. J. Košir, S. Mousavihashemi, B. P. Wilson, E. L. Rautama and T. Kallio, Solid State Ionics, 2022, 380, 115943.
46. A. Dive, K. Kim, S. Kang, L. F. Wan and B. C. Wood, Phys. Chem. Chem. Phys., 2024, 26, 13762–13772.
47. S. Monismith and J. Qu, J. Phys. Chem. C, 2021, 125, 10777–10785.
48. S. Monismith and J. Qu, Extreme Mech. Lett., 2022, 52, 101671.
49. A. Stukowski, Modell. Simul. Mater. Sci. Eng., 2012, 20, 045021.
50. A. Stukowski, Modell. Simul. Mater. Sci. Eng., 2010, 18, 015012.
51. P. M. Piaggi and M. Parrinello, J. Chem. Phys., 2017, 147, 114112.
52. R. K. R. Addula, S. K. Veesam and S. N. Punnathanam, Mol. Simul., 2021, 47, 824–830.
53. J. L. González-Velázquez, in A Practical Approach to Fracture Mechanics, ed. J. L. González-Velázquez, Elsevier, 2021, ch. 3, pp. 75–106,  DOI:10.1016/B978-0-12-823020-6.00003-7.
54. J. L. Allen, J. Wolfenstine, E. Rangasamy and J. Sakamoto, J. Power Sources, 2012, 206, 315–319.
55. T. Thompson, J. Wolfenstine, J. L. Allen, M. Johannes, A. Huq, I. N. David and J. Sakamoto, J. Mater. Chem. A, 2014, 2, 13431–13436.
56. L. J. Miara, W. D. Richards, Y. E. Wang and G. Ceder, Chem. Mater., 2015, 27, 4040–4047.
57. L. Cheng, J. S. Park, H. M. Hou, V. Zorba, G. Y. Chen, T. Richardson, J. Cabana, R. Russo and M. Doeff, J. Mater. Chem. A, 2014, 2, 172–181.
58. C. Ma, Y. Q. Cheng, K. B. Yin, J. Luo, A. Sharafi, J. Sakamoto, J. C. Li, K. L. More, N. J. Dudney and M. F. Chi, Nano Lett., 2016, 16, 7030–7036.
59. C. G. Cole, D. Fincher, C. Roach, R. Osmundsen, B. W. Sheldon, W. Craig Carter, J. LeBeau and Y.-M. Chiang, in Electrochemical Society Meeting Abstracts 245, The Electrochemical Society, Inc., 2024, pp. 2300–2300.
60. C. a. G. Fincher, C. Gilgenbach, R. Osmundsen, C. Roach, A. Penn, M. Thoulass, W. Craig Carter, B. Sheldon, J. Lebeau and Y.-M. Chiang, et al., in Proceedings of 24th International Conference on Solid State Ionics (SSI24), 2024.
61. C. Thomas, W. Zhang, M. R. Chancey, M. Di Michiel, K. Garman, Y. Y. Wang, S. Harris, D. P. Finegan, Y. Q. Wang and C. M. Ban, Cell Rep. Phys. Sci., 2025, 6, 102544.
62. P. Posch, S. Lunghammer, S. Berendts, S. Ganschow, G. J. Redhammer, A. Wilkening, M. Lerch, B. Gadermaier, D. Rettenwander and H. M. R. Wilkening, Energy Storage Mater., 2020, 24, 220–228.
63. R. Dubey, J. Sastre, C. Cancellieri, F. Okur, A. Forster, L. Pompizii, A. Priebe, Y. E. Romanyuk, L. P. H. Jeurgens, M. Kovalenko and K. Kravchyk, Adv. Energy Mater., 2021, 11, 2102086.
64. J. Sastre, A. Priebe, M. Döbeli, J. Michler, A. N. Tiwari and Y. E. Romanyuk, Adv. Mater. Interfaces, 2020, 7, 2000425.
65. J. Sastre, M. H. Futscher, L. Pompizi, A. Aribia, A. Priebe, J. Overbeck, M. Stiefel, A. N. Tiwari and Y. E. Romanyuk, Commun. Mater., 2021, 2, 1–10.
66. K. Kim, A. Dive, A. Grieder, N. Adelstein, S. Kang, L. F. Wan and B. C. Wood, J. Chem. Phys., 2022, 156, 221101.
67. G. Larraz, A. Orera and M. L. Sanjuán, J. Mater. Chem. A, 2013, 1, 11419–11428.
68. F. Tietz, T. Wegener, M. T. Gerhards, M. Giarola and G. Mariotto, Solid State Ionics, 2013, 230, 77–82.
69. M. J. Klenk and W. Lai, Solid State Ionics, 2016, 289, 143–149.
70. A. P. Thompson, H. M. Aktulga, R. Berger, D. S. Bolintineanu, W. M. Brown, P. S. Crozier, P. J. I. Veld, A. Kohlmeyer, S. G. Moore, T. D. Nguyen, R. Shan, M. J. Stevens, J. Tranchida, C. Trott and S. J. Plimpton, Comput. Phys. Commun., 2022, 271, 108171.
71. A. H. Larsen, J. J. Mortensen, J. Blomqvist, I. E. Castelli, R. Christensen, M. Dulak, J. Friis, M. N. Groves, B. Hammer, C. Hargus, E. D. Hermes, P. C. Jennings, P. B. Jensen, J. Kermode, J. R. Kitchin, E. L. Kolsbjerg, J. Kubal, K. Kaasbjerg, S. Lysgaard, J. B. Maronsson, T. Maxson, T. Olsen, L. Pastewka, A. Peterson, C. Rostgaard, J. Schiotz, O. Schütt, M. Strange, K. S. Thygesen, T. Vegge, L. Vilhelmsen, M. Walter, Z. H. Zeng and K. W. Jacobsen, J. Phys.: Condens. Matter, 2017, 29, 273002.
72. R. Freitas, M. Asta and M. De Koning, Comput. Mater. Sci., 2016, 112, 333–341.

---