Atomic-scale mechanisms of interphase formation at lithium–glassy sulfide electrolyte interfaces

Abstract

Glassy sulfide solid electrolytes are promising candidates for all-solid-state lithium batteries owing to their high ionic conductivity and favorable mechanical properties. However, their thermodynamic instability against lithium metal leads to the formation of a complex solid electrolyte interphase (SEI), whose formation mechanisms remain poorly understood. Here, we employ machine-learning force-field molecular dynamics simulations to investigate SEI formation at Li metal interfaces with three representative glassy sulfide electrolytes: 50Li₂S–50SiS₂ (LiSiS), 60Li₂S–32SiS₂–8P₂S₅ (LiSiPS), and 75Li₂S–25P₂S₅ (LiPS). Our simulations reveal that SEI growth follows a power-law dependence across all compositions, with faster growth in P-rich systems. Interfacial reactions proceed through preferential decomposition of P–S and Si–S structural units, with phosphorus exhibiting more rapid reduction kinetics than silicon. The resulting SEI is dominated by an amorphous Li₂S-rich phase, whose composition and transport properties depend strongly on electrolyte chemistry. Notably, a stochastic crystallization event is observed in LiPS, forming a defect-rich, P-containing Li₂S phase that strongly slowed SEI thickening. These findings provide atomic-scale insights into the interplay between glass composition, reaction kinetics, and interphase stability, offering guidance for the rational design of stable lithium–electrolyte interfaces in solid-state batteries.

---

1. Introduction

All-solid-state lithium-ion batteries (ASSBs) are widely regarded as next-generation energy storage systems due to their potential for high energy density and improved safety and stability compared to conventional Li-ion batteries.^1,2 A key component in ASSBs is the solid-state electrolyte (SSE), which transports Li⁺ ions between electrodes and thus largely determines the battery's performance and viability. Among various SSE candidates, sulfide-based solid electrolytes, including LGPS-type, argyrodite-type, and thiophosphate-type materials, have attracted significant attention for their superior ionic conductivity³ (around 10⁻³ to 10⁻² S cm⁻¹ at ambient conditions) and favorable mechanical properties (soft, ductile nature that ensures good interfacial contact).⁴ These advantages, combined with ease of cold-press processing, make sulfide SSEs promising for practical ASSBs. However, a critical challenge is their limited electrochemical stability window: sulfide electrolytes are only thermodynamically stable up to about 1.6–2.3 V vs. Li/Li⁺.^5,6 In other words, when in contact with a lithium metal anode (0 V), these SSEs are intrinsically prone to reduction. This inherent instability leads to spontaneous interfacial reactions that decompose the electrolyte, yielding products such as Li₂S and Li₃P at the Li/SSE interface.^6,7 The accumulation of these decomposition products forms a solid electrolyte interphase (SEI) – a passivation layer that can electronically isolate the electrolyte from the Li metal. The properties of this SEI layer (its composition, morphology, ionic/electronic transport) critically influence cell performance and lifetime, as it can either protect the SSE and enable stable cycling or undergo continual growth and degradation. Indeed, a stable, Li₂S-rich SEI can serve as a beneficial passivation barrier that mitigates further electrolyte decomposition,^8,9 whereas uncontrolled SEI growth or electronically conductive byproducts (e.g. Li–alloy species) may lead to increased interfacial resistance, lithium dendrite propagation, or cell failure. The formation of the SEI can also couple with electro-chemo-mechanical degradation,¹⁰ leading to stress buildup, void and crack formation, accelerated SSE degradation, and dendrite growth.¹¹ Ensuring the formation of a stable SEI at the Li/SSE interface is therefore widely recognized as a key requirement for high-performance and long-lived ASSBs.⁷

Studying SEI formation at solid–solid interfaces, however, poses unique challenges. Unlike the 10–50 nm organic/inorganic layer found in liquid-electrolyte cells,^12,13 the SEI in ASSBs results from direct chemical reactions between two solids. This solid–solid interfacial SEI is often buried between the anode and SSE, making it difficult to observe or characterize in situ. Experimental techniques like X-ray photoelectron spectroscopy (XPS), transmission electron microscopy (TEM), or time-of-flight secondary ion mass spectrometry (TOF-SIMS) can identify final reaction products or interdiffusion profiles,¹⁴ but capturing the dynamic evolution of the SEI remains extremely difficult, especially for the early-stage nucleation and growth. Moreover, sulfide SSEs can be multi-component materials and, in the case of glassy electrolytes,^15,16 lack long-range order – factors which complicate predicting how the interphase develops. These complexities motivate computational approaches to unravel the atomic-scale mechanisms of interfacial decomposition and SEI formation in ASSBs.

Computational studies have indeed emerged as powerful tools to understand Li/SSE interface chemistry.¹⁷ Ab initio calculations, particularly density functional theory (DFT), have provided valuable thermodynamic insights. For example, DFT calculations have shown that lithium thiophosphate electrolytes, such as Li₃PS₄ or argyrodite-type, are thermodynamically unstable in contact with Li metal, and tend to form Li₂S and Li₃P as decomposition products.^7,18 These calculations reveal the driving forces and possible reaction pathways for SSE reduction. Notably, first-principles thermodynamic analysis suggests that the apparent stability of some SSEs observed experimentally is not intrinsic – rather, it arises because the initial decomposition products create a passivating interphase that slows further reaction.⁶ In other words, the SSE's “stable” voltage window is extended by kinetic stabilization: the formation of an SEI (analogous to the liquid-electrolyte SEI) can mitigate the extreme chemical potential difference at the interface and protect the remaining electrolyte. While DFT is useful for mapping out such reactions and interphase energetics, its high computational cost limits simulations to small system sizes (tens to hundreds of atoms) and ultrashort timescales (picoseconds). As a result, purely DFT-based simulations cannot fully capture the realistic morphology or long-time dynamical evolution of SEI. For instance, modeling a realistic Li/SSE interface requires thousands of atoms and nanoseconds of simulation to observe phases nucleating and growing – a regime far beyond standard DFT capabilities.

Classical atomic force field methods offer a technique to simulate systems at larger length and time scales. With parametrized interatomic potentials, one can model millions of atoms over nanoseconds to microseconds, enabling, for example, studies of ion transport across entire cells or long-term structural relaxations. Sun et al. utilized classical force field with constant potential method to investigate the ion transport in argyrodite-type Li₆PS₅Cl constrained by two electrodes under external potential.¹⁹ This allows the study of solid electrolytes involving interfacial reactions with electrodes. However, a fundamental limitation of most classical force fields is their inability to handle arbitrary chemical reactions. The potential functions are typically fitted to a fixed network of bonds or expected coordination environments. Standard force fields lack the flexibility to form new phases or compounds in situ unless those reactions are pre-parametrized (as in specialized reactive force fields like ReaxFF). Even when using reactive force fields, achieving quantitative accuracy for a complex multi-element system (Li–P–S–Si in our case) is extremely challenging, and developing a reliable parameter set requires enormous effort. Thus, conventional molecular dynamics (MD) alone has not been able to fully capture the SEI formation process at Li/SSE interfaces with high fidelity.

Machine learning force fields (ML-FFs) have recently emerged as a promising solution to bridge the accuracy and scale gap in battery interface simulations. ML-FFs (also called neural network potentials) are trained on extensive DFT data, enabling near-DFT accuracy in calculating energies and forces while retaining computational speed comparable to classical MD. Crucially, ML-FF models do not require predefined bonding rules – they learn the potential energy surface directly from electronic structure calculations. This means an ML-FF can naturally accommodate bond breaking and formation events, allowing chemical reactions (like SSE decomposition) to occur spontaneously in MD simulations.²⁰ Several recent works demonstrate the power of ML-FFs in battery interface research. Ren et al. investigated the Li|Li₃PS₄ interfaces formation using ML-FF with model sizes of 40 nm, revealing a four-stage evolution process of the SEI formation beginning with the initial fast ion diffusion stage, followed by nucleation and the growth of Li₂S, and stabilization.²¹ Chaney et al. examined the SEI formation of the Li|Li₆PS₅Cl system, proposing a two-step growth mechanism: the formation of an amorphous phase, followed by the kinetically slower crystallization of the reaction products into a 5Li₂S·Li₃P·LiCl solid solution.²² These simulations, achieved with accuracies near that of DFT, provided unprecedented atomic-scale insight into how SEIs develop on crystalline sulfide electrolytes.

Despite these advances, comparatively less attention has been paid to glassy solid electrolytes (GSEs) in the context of Li metal interfaces. Most computational investigations so far have focused on crystalline or ordered SSEs, yet many practical sulfide electrolytes are amorphous or glass-ceramics. For example, the Li₂S–P₂S₅ glass family is a commercially important SSE. Glassy SSEs differ fundamentally from crystalline materials in their structure and potentially in their interfacial chemistry. The lack of long-range order in a glass indicates the Li penetration and reaction may not follow the same pathways observed in crystals in which specific facets or grain boundaries dictate reaction fronts.²³ Moreover, glassy electrolytes often contain mixed anion/cation frameworks – for instance, adding SiS₂ into Li₂S–P₂S₅ forms a ternary glass with both P and Si as network formers.²⁴ The presence of different glass-forming units could lead to preferential decomposition of one unit over the other, altering the SEI composition and its growth kinetics. For example, P–S bonds might break more readily than the Si–S bonds under reducing conditions, as Si-rich GSE compositions show higher weighted average bond dissociate enthalpy,²⁵ yielding different rates of Li₂S formation. Indeed, recent work in Li metal batteries has suggested that amorphous electrolytes and crystalline electrolytes exhibit different interfacial reaction dynamics.^12,26 However, a detailed mechanistic understanding of SEI formation at Li/glassy-SSE interfaces is still lacking, making this an important open question for the field. Addressing this gap is crucial, since glassy sulfide electrolytes (such as the Li₂S–P₂S₅ system) are among the most attractive SSEs for practical ASSBs – their high ionic conductivity and lack of grain boundaries (which helps prevent dendrite penetration) are significant advantages,^15,24 if their interface with lithium can be made stable.

In the present work, we develop an ML-FF for describing Li–Si–P–S interfacial problems and apply it to investigate the interfacial chemistry between lithium metal and sulfide-based glassy solid electrolytes. We focus on the prototypical Li–Si–P–S glass systems with three representative compositions: a Si-rich binary glass 50Li₂S·50SiS₂ (LiSiS), a mixed P + Si ternary glass 60Li₂S·32SiS₂·8P₂S₅ (LiSiPS), and a P-rich binary glass 75Li₂S·25P₂S₅ (LiPS). The numbers preceding each compound denote mole percentage. These compositions allow us to examine the effects of network formers (SiS₂ vs. P₂S₅) and glass structural disorder on SEI formation. Using large-scale ML-FF based MD simulations at 300 K, we directly simulate the initial contact between Li metal and each glassy SSE and track the evolution of the interphase over tens of nanoseconds. By comparing the three systems, we show how compositional complexity and the absence of long-range order impact the SEI growth mechanism and stability. Our simulations reveal that all three glasses undergo an interfacial reaction whereby the PS₄ (phosphorus-sulfide) and SiS₄ (silicon-sulfide) structural units in the glass are reduced, forming an amorphous Li₂S-rich layer at the interface. In addition, the SEI growth is diffusion-limited, but its rate and morphology differ by compositions. Notably, P-containing glasses (LiSiPS and LiPS) show faster SEI growth kinetics than the Si-only glass (LiSiS), which we correlate with the more facile breakage of P–S bonds compared to Si–S bonds in the alloyed network. In the P-rich LiPS composition, the SEI's Li₂S-like domains become sufficiently developed that a partial crystallization to a Li₂S phase is observed within our MD simulation time (nanoseconds) – a phenomenon not observed in the Si-containing glass systems. These findings provide new atomic-scale insights into how structural disorder and glass chemistry influence interphase formation in sulfide SSEs.

2. Computational methodology

2.1 Development of the machine learning force field

The ML-FF was developed using the DeepPot-SE framework,^27–29 adopting the same hyperparameters as in our previous study, including a cutoff radius of 6.0 Å for local atomic environments.³⁰ The architecture consists of a three-layer embedding network (25, 50, 100 neurons) and a three-layer fitting net with 120 neurons per layer. To ensure that the model accurately captures interfacial environments and reaction behaviors, the training dataset was significantly expanded relative to our previous work.³⁰ The additional training data includes lithium surface structures (bcc (100, 110,111), fcc (100, 111), hcp (0001) orientations with >10 Å vacuum layer), binary LiP and LiSi bulk structures (Li₃P, Li₃P₇, LiP, LiP₅, LiP₇, Li₂Si, Li₇Si₂, and LiSi), and Li|SSE interface structures in which we focus on Li bcc(001) interface with amorphous SSE. The interface structures were generated using pymatgen³¹ and dpdata³² packages.

The training dataset was iteratively enriched through 64 iterations of active learning using the DP-GEN framework.³³ In each iteration, configurations were sampled from both canonical (NVT) and isothermal–isobaric (NPT) MD simulations over a temperature range between 300 and 1000 K. Then the DFT calculations were performed on the selected configurations using the Vienna Ab initio Simulation Package (VASP).^34–36 The PBEsol exchange-correlation functional,³⁷ a plane-wave energy cutoff of 500 eV, and a Monkhorst–pack k-point spacing of 0.5 Å⁻¹, were employed in DFT simulations, which are consistent with our previous work.³⁰ The final training dataset contains 215 095 configurations with corresponding total energies, atomic forces, and virial stress tensors. With the final training dataset, we carried out 8 million training steps, during which the learning rate was gradually reduced from 2 × 10⁻³ to 3.5 × 10⁻⁸.

2.2 Molecular dynamics simulations and analysis

MD simulations were performed using the LAMMPS package³⁸ with the developed ML-FF. The Li|SSE interface models were constructed by placing a Li metal slab in contact with a glassy SSE slab. The Li slab was cleaved along the bcc (001) direction and then was relaxed, consisting of 3072 atoms with a thickness of 83 Å. The glassy SSE models were generated via the melt-quench method: initial atomic configurations were generated by randomly putting molecular fragments (Li, PS₄, P₂S₆, P₂S₇, SiS₄, Si₂S₆, Si₂S₇, Si₄S₁₀) in the supercell via PACKMOL³⁹ package, and then equilibrated at 1000 K using isothermal–isobaric (NPT) ensemble to obtain fully molten structures. Next, the molten SSEs were quenched to 300 K within 500 ps to obtain the glass phase and then they were deformed to match the lateral dimensions (x and y) of the Li slab. Finally, the glass structures were equilibrated at 300 K for 10 ps using canonical (NVT) ensemble before contacting the Li slab. Three representative compositions were investigated: 50Li₂S–50SiS₂, 60Li₂S–32SiS₂–8P₂S₅, and 75Li₂S–25P₂S₅, hereafter denoted as LiSiS, LiSiPS, and LiPS, respectively. To account for structural variability inherent to amorphous systems, six statistically independent SSE structures were generated for each composition using different random seeds for the initial velocity assignment in the melt-quench process.

Periodic boundary conditions (PBCs) were applied in all three spatial directions. To restrict the dynamics to a single effective interface and prevent spurious interactions across periodic images, the top and bottom 7 Å atoms were ‘frozen’ along the z axis during MD simulations and excluded from all subsequent analyses. MD simulations on the Li|SSE interface models were conducted in the NPT ensemble at 300 K with a timestep of 1.0 fs. For each of the 18 systems (3 different compositions × 6 amorphous systems), 15 ns simulations were carried out.

To quantitatively determine the interface positions and characterize the SEI region, we fitted logistic functions to atomic density and atomic potential energy profiles along the surface normal (z) direction. For interfacial profiles with one step, a standard logistic function was used:

For double-step profiles, a sum of two logistic functions was employed:

where y_i, x_i, w_i are fitting parameters. The interface transition widths were defined as the region over which property changes from 10% to 90% of the full step height (y₂ − y₁ or y₃ − y₂), and the SEI thickness was calculated as the distance between two transition centers x₂ − x₁.

Coordination number (CN) profiles were computed by binning atoms along the z-direction. For each atom pair, neighbors were counted within element-specific cutoffs corresponding to the first minimum of the partial radial distribution functions (RDFs): P–Li (2.86 Å), Si–Li (3.49 Å), S–Li (3.39 Å), P–S (2.73 Å), Si–S (2.80 Å), Li–Li (4.00 Å), and Si–Si (2.64 Å). Density profiles were calculated analogously using atomic mass and bin volumes along the z-axis. The self-diffusion coefficient D was calculated from the mean squared displacement (MSD) of selected atoms using , where d is the dimensionality (d = 1 for diffusion along the z-direction, and d = 2 for diffusion in the xy-plane), N is the number of atoms considered, and t is the elapsed time. Atoms belonging to the fixed boundary layers were excluded from the diffusion analyses, as well as CN and density analyses.

3. Results and discussion

3.1 Validation of ML-FF

Fig. 1 shows the performance of the developed ML-FF used in this study compared to DFT reference data. The ML-FF achieves a root mean squared error (RMSE) of 11 meV per atom for energies and 0.2 eV Å⁻¹ for forces, showing slightly degraded performance compared to our previous ML-FF trained only on bulk datasets. This is attributed to the broader distribution of energies and energy gradients encompassed in the new training data. Nevertheless, the ML-FF still maintains good overall accuracy. In addition, we compared the Li|LiSiPS interfacial structures obtained from the AIMD and ML-FF trajectories. The good agreement between the RDFs (Fig. S1 in the SI) further supports the validity of the developed ML-FF.

Fig. 1 Comparison between the ML-FF and DFT reference data. The reference dataset includes configurations selected from 3000 structures spanning 305 distinct systems. (a) Total energies and (b–d) atomic force components along the x, y, and z directions (f_x, f_y, and f_z, respectively).

3.2 SEI formation kinetics and growth mechanism

Fig. 2 displays the simulation snapshot for the 50Li₂S–50SiS₂, 60Li₂S–32SiS₂–8P₂S₅, 75Li₂S–25P₂S₅ compositions at 1 ns, along with their density and atomic potential energy (PE) profile. Three distinct regions are observed on the snapshots (Fig. 2a–c): the Li metal region at the bottom, the solid electrolyte interphase region at the middle, and the unreacted solid-state electrolyte region at the top. Fig. 2d–f presents the atomic PE of Li (E^Li_i) along the z-direction. E^Li_i in metallic Li is −1.90 eV, which is higher than the corresponding values in the SEI (−2.46 to −2.41 eV) and the glass electrolyte (−2.78 to −2.75 eV), indicating a stronger local energetic instability as Li contacts electrolyte. From a thermodynamic perspective, the E^Li_i can be viewed as a proxy for the Li chemical potential, reflecting the driving force of forming the SEI between Li and amorphous electrolyte. On the other hand, E^Li_i is a scalar descriptor encoding the local environment of the surrounding atoms. Fig. 2g–i shows the density profiles of the three compositions, revealing a pronounced change in density at the formed SEI layer (∼z = 70 Å). In contrast to the E^Li_i profile, which exhibits an additional discontinuity at the SEI|SSE interface, the density profiles display only a single distinct jump at the Li|SEI interface across all three systems. This indicates that the energetic heterogeneity across the SEI|SSE boundary is not accompanied by a corresponding density contrast. The SEI formation primarily alters the local packing density near the Li surface, while the SEI|SSE interface remains structurally continuous.

Fig. 2 Snapshots of simulation trajectories and interfacial characteristics. (a–c) Representative snapshots at 0 ns and 1 ns for LiSiS (a), LiSiPS (b), and LiPS (c), illustrating interfacial evolution during simulation. (d–f) Atomic potential energy profiles of Li along the z direction for the three systems. (g–i) Corresponding density profiles along the z direction, highlighting structural variations across the Li|SEI|SSE interfaces. The Li, Si, P, and S atoms are represented by purple, beige, orange, and brown spheres, respectively.

Using E^Li_i profiles in Fig. 2(d–f), we quantify the interfacial thicknesses of both the Li|SEI and SEI|SSE regions. At 1.0 ns, the Li|SEI interface in the LiPS SSE exhibits the narrowest transition region (∼4.28 Å), whereas the LiSiS and LiSiPS SSEs display slightly broader transition widths of 5.09 Å and 5.62 Å, respectively. In contrast, the SEI|SSE interface is substantially more diffuse. LiPS shows a relatively narrow transition region of ∼10.5 Å, while LiSiS and LiSiPS exhibit much broader, smeared interfacial profiles with thicknesses of ∼22 Å. Notably, the thicknesses of both the Li|SEI and SEI|SSE interfaces remain largely unchanged beyond 1 ns.

Next, we examine the time evolution of the SEI thickness for the three SSE compounds (LiSiS, LiSiPS, LiPS). By identifying the midpoint positions of the Li|SEI and SEI|SSE interfaces, we determine the SEI thickness as the distance between these two boundaries during the MD simulations. Fig. 3 presents the SEI thickness evolution over the first 1.0 ns, while Fig. S2 in the SI shows the corresponding results over the initial 5 ns. All three glassy SSE compositions exhibit an initial spontaneous formation of an amorphous interphase (SEI) when contacted with Li metal. The SEI growth is rapid at early times and then slows, following a power law time-dependence model. This behavior suggests that Li⁺ transport into the SSE (and across the formed SEI) governs the reaction rate, as freshly reduced phases accumulate at the interface and progressively impede further Li diffusion. Indeed, we find that Li metal infiltrates into the glassy SSE and drives reduction reactions at a moving reaction front, analogous to the initial “fast ion diffusion” stage reported in previous ML-FF studies of crystalline SSE interfaces.²¹ Quantitatively, the SEI thickens fastest in the P₂S₅-rich glass (LiPS) and slowest in the SiS₂-only glass (LiSiS), with the mixed composition (LiSiPS) intermediate – i.e. the growth rates decrease in the order LiPS > LiSiPS > LiSiS. This trend indicates that composition strongly influences reaction kinetics and is consistent with different reduction behavior of P–S and Si–S units, as discussed below. Nevertheless, in all cases the power-law behavior persists, suggesting broadly similar transport-limited SEI growth behavior across these glasses.

Fig. 3 Evolution of SEI thickness during initial growth. SEI thickness as a function of time for the three SSE compounds (LiSiS, LiSiPS, and LiPS) during the initial growth stage up to 1 ns. The shaded region represents the uncertainty estimated from six independent MD simulations.

The diffusion-limited growth is consistent with prior theoretical and experimental reports that an insulating SEI layer self-limits its own growth by blocking electron transfer while still allowing Li⁺ transport.⁴⁰ Notably, first-principles calculations have shown that sulfide SSEs (e.g. Li₃PS₄) will thermodynamically decompose to Li₂S and Li₃P at a Li interface, but that these products are electronically insulating and thus passivate further reaction (kinetically stabilizing the interface).⁴⁰ Our simulations confirm this picture: an initial burst of interfacial reactions produces a Li₂S-rich SEI that slows subsequent decomposition. Such behavior was also observed by Wenzel et al. in experiments, where a thin reduced layer (Li₂S–Li₃P) formed at a Li/Li₇P₃S₁₁ interface and was found to be low in interfacial resistance and remarkably stable.⁴¹

3.3 Reaction pathways and interphase composition

To gain mechanistic insight into the reactions underlying SEI formation, we analyzed the radial distribution functions (RDFs) for representative atomic pairs in the SEI and SSE regions (Fig. S3 and S4 of SI) and performed time-dependent coordination analyses of P–S and Si–S bonds (Fig. 4). At the atomic level, the SEI composition evolves via the breakdown of the SSE's PS₄ and SiS₄ structural units (in P-containing and Si-containing glasses, respectively). As shown in Fig. S2 of SI, the P–S peak shows a clear shift from ∼2 Å in SSE to ∼4 Å in SEI for both LiPS (Fig. S3d) and LiSiPS (Fig. S3b) simulations, suggesting the breakage of P–S bonds as SEI forms. The intensity of the first Si–S RDF peak also shows a small value in the SEI compared to SSE, while a much higher peak is present between 4.2–4.5 Å in SSE for both LiSiS and LiSiPS compounds, indicating the decomposition of SiS₄ tetrahedra. Furthermore, the coordination numbers of P–S and Si–S pairs exhibit different reaction dynamics. The coordination number and RDF results show more extensive P–S bond breakage than Si–S bond cleavage under these conditions, as shown in Fig. 4 and S3. Within only a few nanoseconds, essentially all P–S tetrahedra in the SEI region have been reduced: sulfur is stripped from P₂S₅ units to form S²⁻ anions (ultimately Li₂S), and P is freed from the framework to form Li₃P. By contrast, the Si–S bonds prove more robust – in the thiosilicate glass (LiSiS), only those SiS₄ units in the immediate vicinity of Li metal become fully decomposed during the same timeframe. Thus, the P-containing glasses undergo more complete and rapid anion reduction, whereas the Si-rich glass retains partially intact SiS₄ units just beyond the SEI|SSE interface. This trend is consistent with composition-dependent kinetics in Fig. 3: P-containing units appear to provide faster reduction pathways, whereas Si-containing units are more persistent. It also means that the resulting SEI compositions differ. In all cases, Li₂S is the dominant inorganic phase in the SEI (as expected from S²⁻ formation).

Fig. 4 Coordination number profile along z-axis as a function of time for simulations of Li|LiSiPS (60Li₂S–32SiS₂–8P₂S₅).

As shown in Fig. S4, the first RDF peak of Li–Li pair shifts from 3.5 Å in the SSE region to 3.0 Å in the SEI region. This shortened distance corresponds to the Li–S–Li distance in Li₂S. In contrast, the Li–S pair in RDFs remains relatively unchanged. This together indicates the formation of Li₂S-like short-range ordering within the SEI layer. In contrast, the network-forming elements behave differently. Phosphorus is mostly converted to reduced P species (such as Li–P compounds) across the SEI region, while silicon remains behind largely as unreacted or partially reduced SiS_x units near the SSE-side. The Li–Si pair coordination number shows a high value near the Li|SEI interface, in contrast to the near uniformed result of P–S pair across the whole SEI region (Fig. S5–S7 of SI). This suggests that the formation of Li–Si alloys from SiS₂ is kinetically hindered at room temperature, in stark contrast to analogous crystalline SSEs like LGPS (Li₁₀GeP₂S₁₂) where in situ XPS and DFT studies have shown Li readily reduces Ge(IV) to Li–Ge alloys alongside Li₂S and Li₃P.⁴²

The compositions of the SEI layers formed on different SSE compounds (Table 1) remain relatively stable throughout the simulations. For LiPS, the SEI composition is approximately Li : P : S = 67% : 7% : 26%, which is close to the theoretical thermodynamic decomposition ratio of 68.75% : 6.25% : 25% predicted from the reaction Li₃PS₄ + 8Li → 4Li₂S + Li₃P. In contrast, for the LiSiS and LiSiPS solid electrolytes, the SEI compositions exhibit excess sulfur and lithium deficiency relative to the expected decomposition products (Li₂₁Si₅–Li₃P–Li₂S). This deviation indicates that the Li_xSi phase in the SEI is less lithiated (Li-poor) and contains partially reduced Si–S units.

Table 1 Elemental number densities of Li, Si, P, and S in the SEI region for the three SSE compounds at 1 ns of simulation time. The atoms are counted between the Li|SEI and SEI|SSE interfaces and normalized by the SEI region volume

Composition	Li (Å^−3)	Si (Å^−3)	P (Å^−3)	S (Å^−3)
50Li2S–50SiS2	0.0391	0.0052		0.0169
60Li2S–32SiS2–8P2S5	0.0375	0.0033	0.0017	0.0178
75Li2S–25P2S5	0.0419		0.0041	0.0166

---

Overall, our findings are consistent with the overall picture of the sulfide reduction mechanisms (formation of Li₂S, Li₃P, and Li_xSi species) that have been predicted by previous first-principles studies and observed experimentally.⁴⁰ At the same time, this study reveals that glassy vs. crystalline SSEs can follow distinct reaction pathways: the disordered network of a glass allows a more homogeneous, gradual decomposition (with some bonds surviving deeper in the glass), whereas crystalline SSEs often undergo abrupt decomposition along specific lattice planes or defect sites.⁸ Indeed, interfaces with crystalline β-Li₃PS₄ are known to generate structural defects and heterogeneities upon reduction,⁸ whereas the amorphous structure of the glassy SSE yields a more continuous SEI layer without the lattice-mismatch dislocations inherent to crystals. This intrinsic difference likely leads to the more uniform SEI growth in glasses, in contrast to the facet-dependent or grain-boundary localized reactions that can occur in polycrystalline sulfide electrolytes.

3.4 Crystallization and structural evolution of the SEI

A notable finding from our simulations is that the initially amorphous SEI can undergo partial crystallization under ambient conditions, although this process is highly sensitive to composition and local structural fluctuations. In five out of six independent simulations for LiPS composition, the SEI remained amorphous over the 15 ns simulation. However, in one independent MD run, a spontaneous nucleation event occurred around ∼4 ns, leading to the formation of a crystalline Li₂S-rich phase within the SEI. This crystallization was evidenced by the emergence of long-range order in the S sublattice and a plateau in SEI thickness, as shown in Fig. 5(a–c). In contrast to LiPS, all six independent MD simulations for the other two SSE compounds (LiSiS and LiSiPS) exhibit exclusively amorphous SEIs over the 15 ns simulation period.

Fig. 5 Snapshots and evolution of the crystallized SEI structure. (a) Representative snapshot of the LiPS compound system at 5 ns, showing the formation of a crystallized SEI. (b) Magnified view of the crystallized SEI region outlined in red in (a). (c) Atomic positions in the crystallized SEI region, superimposed into the cubic cell outlined in green in (b). (d) Time evolution of SEI thickness over 10 ns for all six independent Li|LiPS independent simulations. Run 1 exhibits crystalline SEI, whereas the other five runs remain amorphous. The Li, P, and S atoms are represented by purple, orange, and brown spheres, respectively.

In the simulation exhibiting SEI crystallization, the crystalline domain consists predominantly of Li₂S with measurable incorporation of phosphorus into the lattice. Part of the phosphorus released from P₂S₅ is incorporated as dopants or substitutional species within the Li₂S framework (Fig. 5b and c), rather than segregating into a separate Li₃P phase. Consequently, no crystalline Li₃P or other P-rich phases are observed; instead, phosphorus remains either amorphous or in solid solution within the Li₂S matrix. Once this Li₂S-rich crystalline domain formed in the SEI, the thickness of SEI (Fig. 5d) showed only a negligible increase over the simulation time window (∼0.0045 Å ns⁻¹). This suggests that, in this simulation and over the simulated time range, the crystalline Li₂S-rich region may partially passivate the interface and slow further chemical growth. Interestingly, the Li|SEI boundary still migrated slightly toward the Li anode (at ∼0.04 Å ns⁻¹) even after SEI crystallization (Fig. S8), implying ongoing slow Li penetration via defects but no net thickening of the SEI.

To further test reproducibility for this SEI crystallization event, we performed five additional 10 ns simulations from the exact same initial LiPS structure but using different random seeds for initial velocity creation. In these simulations, two out of five trajectories also exhibit clear SEI crystallization (Fig. S9), although the nucleation location and time differ. Taken together, these results indicate that crystallization is a stochastic pathway that can recur in some local environments. However, due to the limited simulation cell size (∼3 × 3 × 18 nm³) and time window (<15 ns), they are not sufficient to establish as the general or dominant behavior of sulfide glassy solid electrolytes. The stochastic nature of the crystallization – occurring in only one out of several runs within MD timescale – points to the role of chemical and structural fluctuations in nucleating an ordered phase within an amorphous matrix. In a glassy electrolyte, the SEI may need a favorable cluster or composition (e.g. a sufficiently large Li₂S-rich region) to trigger crystallization, which was realized in that particular LiPS MD run. The relatively high P content in LiPS may facilitate this process by generating a larger fraction of Li₂S through complete reduction of P₂S₅ and by promoting faster ionic transport that enables structural rearrangement. In contrast, the Si-containing systems (LiSiPS and LiSiS) do not crystallize within the MD simulation timescale. Their SEIs remain disordered within the simulated time range, which might be related to the persistent presence of partially reduced Si–S units and a lower volume fraction of pure Li₂S.

These results align with previous simulations on crystalline SSEs, where a two-step SEI formation was observed: an initial amorphous layer formation followed by a slower crystallization of decomposition products. Chaney et al., for example, reported that the Li/Li₆PS₅Cl interface first forms an amorphous interphase, and upon extended simulation it crystallizes into a mixed 5Li₂S·Li₃P·LiCl phase (a solid solution).²² Our findings are consistent with this two-step mechanism; however, in the glassy systems the crystallization step is not as deterministic. When it does occur (as in LiPS), the end product is essentially Li₂S with dopants rather than a well-defined multi-component crystal. This could be viewed as an anion-disordered variant of the sulfide SEI – in LiPS, lacking a third component like Cl, the system crystallized into a Li₂S-based phase with P substituents, whereas in Li₆PS₅Cl the presence of Cl leads to a ternary solid solution crystallite. Importantly, the onset of crystallization in our LiPS simulation coincides with a transition to a stabilized SEI thickness, reinforcing that long-range ordering of Li₂S can act to arrest further reaction (by creating an electron-insulating, ion-saturated barrier). This finding supports previous proposals that thin Li₂S interlayers can passivate the Li/SSE interface.

From a practical perspective, our results suggest that, under favorable conditions, a glassy SSE interface may evolve toward a more stable, Li₂S-rich crystalline interphase. However, the persistence of amorphous SEIs in most simulations indicates that crystallization in glasses likely requires extended time scales or external stimuli, such as elevated temperature or applied bias or stress. More robust quantification of crystallization probability and long-time passivation behavior will require longer simulations, larger interfacial cells, and more quantitative structural analysis, potentially combined with enhanced-sampling methods. Finally, while partial crystallization may enhance mechanical integrity and chemical stability, it could also introduce interfacial stresses or brittleness associated with volume changes. Therefore, crystallization of glass-derived SEIs may represent a double-edged sword for long-term interfacial stability in solid-state batteries.

3.5 Ion mobility and defect chemistry in SEI region

Despite being composed of electronically insulating compounds, the SEI layers formed here remain conducive to Li⁺ transport – a critical requirement for stable cycling. Fig. 6 shows the mean-squared displacement (MSD) of Li, Si, P, and S atoms along the in-plane (XY, parallel to the interface) and out-of-plane (z, perpendicular to the interface) directions. All four species exhibit diffusive motion, with overall displacement dominated by the Z direction. For Li ions, diffusion is only weakly anisotropic: the diffusion coefficient along Z (D_z) is comparable to that in the XY plane (D_xy), being only about twofold higher. Specifically, the D_z/D_xy ratios are 2.34, 2.28, and 2.39 for the LiSiS, LiSiPS, and LiPS systems, respectively, as summarized in Table 2.

Fig. 6 Mean squared displacement (MSD) of Li, Si, P, and S atoms in three SSE compounds. (a–d) MSD of Li, Si, P, and S atoms in the lateral directions (X and Y) parallel to the interface plane. (e–h) MSD of Li, Si, P, and S atoms in the vertical direction (Z) perpendicular to the interface.

Table 2 The calculated diffusion coefficients for Li, Si, P, and S species in lateral directions (X or Y) parallel to the interface plane and vertical direction along the interface normal (Z direction) in three SSE compounds from MD simulations

Atoms	Composition	Dxy (10^−7 cm^2 s^−1)	Dz (10^−7 cm^2 s^−1)	Dz/Dxy
Li	50Li2S–50SiS2	1.79	4.16	2.34
	60Li2S–32SiS2–8P2S5	2.30	5.25	2.28
	75Li2S–25P2S5	2.30	5.51	2.39
Si	50Li2S–50SiS2	0.25	5.23	21.3
	60Li2S–32SiS2–8P2S5	0.27	6.79	25.0
P	60Li2S–32SiS2–8P2S5	0.36	7.71	21.2
	75Li2S–25P2S5	0.27	9.92	37.2
S	50Li2S–50SiS2	0.15	4.60	30.4
	60Li2S–32SiS2–8P2S5	0.21	6.43	30.5
	75Li2S–25P2S5	0.25	9.94	39.5

---

Once the SEI in LiPS crystallized into a Li₂S-rich phase, we calculated its Li⁺ ionic conductivity. The Li⁺ in the crystalline SEI shows an average ionic conductivity of 0.397 mS cm⁻¹ (with XY direction 0.329 mS cm⁻¹, 0.534 mS cm⁻¹ along Z direction). The high ionic conductivity of the crystalline SEI is consistent with a defect-chemistry interpretation in which P incorporation and associated disorder may introduce lithium vacancies or vacancy-like pathways. Essentially, P doping turns Li₂S – normally a poor ionic conductor – into a defect-rich superionic conductor. Szczuka et al. showed that Li₃P can form solid solution metastable phase with Li₂S with the anti-fluorite structure that are close to the thermodynamic ground state.⁴³ The formed metastable phase shows a higher room-temperature experimental conductivity around 10⁻⁵ to 10⁻⁴ S cm⁻¹.⁴³ In addition, ReaxFF MD simulations by D'Amore et al. support this mechanism, showing that stoichiometric Li₂S has an extremely low room-temperature Li⁺ conductivity (∼5 × 10⁻⁸ S cm⁻¹), but introducing just a few percent of Li vacancies (e.g. 3% vacant sites, which could result from aliovalent P substitution) raises the conductivity by several orders of magnitude (to ∼2.8 × 10⁻⁴ S cm⁻¹).⁹ As shown in Fig. 5c, the P atoms mainly substitute for S atoms at the 4a sites in the Li₂S lattice and show an occupancy of around 18%. Our simulated SEI exhibits a similar trend: with P-induced vacancies (Fig. 5b), the Li₂S-rich interphase behaves as a mixed sulfide-ion conductor, facilitating Li⁺ migration across the interface.

Notably, even in the amorphous SEI, we expect substantial Li⁺ mobility due to the disordered structure – amorphous Li₂S and Li₃P networks inherently contain coordination defects and pathways for ion hopping. In all simulations, Li atoms were observed diffusing through the developing SEI layer, indicating that the interphase does not block Li⁺ transport. This aligns with experimental impedance studies that report low interfacial resistance after an initial SEI is formed.⁴⁰ In essence, the spontaneously formed SEI acts as a self-regulating membrane. It is ionically conductive (permissive to Li⁺) but electronically insulating, thereby allowing continued Li⁺ transfer while preventing further uncontrolled electron-driven reduction of the electrolyte. Such behavior is analogous to the SEI in liquid-electrolyte batteries, though here the SEI is inorganic (Li₂S–Li₃P-based) and forms at a metal/solid interface. We also examined the diffusion pathways of various species (Si, P, S) within the SEI. As shown in Fig. 6 and Table 2, although the displacements of Si, P, and S atoms along the Z direction are comparable in magnitude to that of Li, their in-plane (XY) mobility is substantially lower. For these species, diffusion in the XY plane is approximately 20–40 times slower than along the Z direction. Notably, their mobilities remain higher than their respective self-diffusion coefficients in the bulk SSE. As expected, Li⁺ cations percolate via the Li₂S-rich matrix, finding “highways” through interconnected S²⁻ networks and vacancy clusters. This is consistent with other reports that Li₂S can sustain facile Li⁺ diffusion along certain crystallographic directions or through grain boundaries.⁸

In contrast, the anionic framework elements (S, P, Si) are largely immobilized in the XY plane of SEI – they rearrange locally to form new phases (e.g. sulfide chains, polyphosphides), but we do not observe long-range diffusion of S, P, and Si away from the interface. This is analogous to the bilayer SEI architectures recently reported in crystalline systems – e.g., a ML-FF study of a Li_xSi |Li₆PS₅Cl interface found a Li₂S/LiCl-rich layer in electrolyte-derived region and a Li₂S/Li₃P-rich layer in anode-derived region, attributed to concentration gradients during SEI growth.⁴⁴ In our case, the Si-containing composition exhibit larger element concentration gradient across the SEI region compared to P-only composition (LiPS). We note that the glassy nature of the SSE yields an SEI that is itself largely amorphous (unless crystallized), meaning the Li⁺ diffusion does not rely on crystalline pathways but rather on the intrinsic free volume and defects of the disordered structure. While the high ionic conductivity of the amorphous SEI lowers interfacial impedance, it also promotes a higher SEI growth rate. Overall, our results illustrate that the interphase formed at Li/glassy–sulfide interfaces is ionically percolating regardless of crystallization or not.

3.6 Comparison to crystalline interfaces and literature

The findings in this work provide a picture of how glassy vs. crystalline SSE interfaces behave differently. In both cases, the thermodynamic outcome is the formation of a Li₂S-based SEI containing reduced P (and other dopants), confirming a common paradigm for Li metal reacting with sulfide electrolytes.⁴⁰ However, the kinetics and structure of the SEI differ substantially between disordered and ordered electrolytes. Crystalline SSEs (such as Li₃PS₄ or Li₆PS₅Cl) tend to decompose at specific crystal interfaces or defect sites, often leading to heterogeneous SEI layers or even crystalline interphases at relatively short times.⁸ For example, simulations of Li|β-Li₃PS₄ have reported a four-stage process with rapid lithium ingress and swift nucleation of Li₂S crystals at the interface.²¹ In Li|Li₆PS₅Cl, the eventual crystallization of a single-phase Li₂S·Li₃P·LiCl interphase was observed, reflecting the long-range ordering tendencies of a system that initially had long-range order.

By contrast, in the glassy SSE interfaces we studied, the SEI initially grows in a more uniform, amorphous manner, without the immediate appearance of distinct crystalline domains. Only one independent MD run of LiPS composition showed a full crystallization event. The structural disorder of the glass may promote a more isotropic reaction front (no crystallographic preference for Li penetration), and the lack of grain boundaries means fewer stress concentrations or fracture pathways during SEI formation. This could be advantageous. We did not observe any void formation or interfacial cracking in our MD simulations, whereas crystalline interfaces often face issues of interfacial voids or mechanical decohesion due to volume changes. In addition, the growth of the amorphous SEI follows a power-law-type behavior. Several time-dependence models have been either observed experimentally or derived theoretically, including t^1/2, t^1/3, tⁿ, and ln t.⁴⁵ In particular, the Wagner-type diffusion model predicts a parabolic growth law d ∝ t^1/2, which has been shown to agree well with experimental observations,⁴⁶ especially in storage/aging experiments.⁴⁵ In our case, the growth of the amorphous SEI region can be well fitted by all four models (R² > 0.95, Table S1 of SI), due to the limited simulation (<15 ns). However, from a predictive perspective, the cubic-root relationship (t^1/3) provides the best overall performance for our systems, as shown in Table S2 of SI.

Experiments have shown that adding a compliant interlayer (like Au) at a Li/crystalline-Li₃PS₄ interface can alleviate voids and improve Li plating uniformity;⁴⁰ in a glassy electrolyte, the intrinsic amorphous structure may similarly accommodate volume changes more gracefully, potentially yielding a more robust contact. On the other hand, a fully amorphous interphase might have lower elastic modulus and be more prone to creep or penetration by Li metal under pressure, whereas a crystalline Li₂S-rich layer (as occasionally formed in LiPS) could provide a harder, though more brittle, barrier. From a chemical stability perspective, both glassy and crystalline sulfide-based SSEs ultimately rely on the formation of a Li₂S-based interphase to attain compatibility with Li metal. Our work extends prior studies on crystalline sulfides by demonstrating that this holds true for multicomponent glassy sulfides as well. For the LiPS GSE, the formed Li₂S–Li₃P interphase exhibit high ionic conductivity while remaining electronically insulating,⁴³ which are ideal for a stable passivating SEI. In contrast, for Si-containing GSEs, the SEI comprises partially reduced Si–S units and partially lithiated Li_xSi phases. The resulting mixed ionic–electronic conducting interphase is therefore intrinsically unstable over long time scales, consistent with experimental observations by Zhao et al.¹⁵ Although these reactive Si environments can temporarily suppress lithium dendrite growth through the reactions Si⁴⁺ + xLi⁺ + xe⁻ → Li_xSi and Li_xSi + yLi⁺ + ye⁻ → Li_x+ySi, continued SEI growth ultimately occurs. This sustained interphase growth arises from the electronic conductivity of the Li_xSi phases, which facilitates ongoing electrochemical reactions at the interface.

These insights provide a mechanistic framework for interpreting experimental observations at Li/anode interfaces involving Si-containing electrolytes or additives. Although direct experimental characterization of SEIs at Li|glass interfaces remains challenging, our results are consistent with trends reported in related systems. For example, Li–P–S glasses, such as 75Li₂S–25P₂S₅, are known to exhibit moderately improved stability against Li metal compared with crystalline thiophosphates, often enabling several stable cycles before failure. This behavior has been attributed to the rapid formation of a passivating Li₂S/P-rich interphase. Our simulations explicitly capture the formation of this passivating layer and quantitatively characterize its thickness, composition, and ionic conductivity. Furthermore, by comparing thiosilicate and thiophosphate glasses, we predict that incorporation of P₂S₅ into sulfide glasses promotes the formation of a more ion-conductive and readily passivating SEI, owing to P-induced structural and electronic defects. In contrast, pure thiosilicate glasses tend to form less conductive interphases dominated by undoped Li₂S, which may lead to higher interfacial resistance. These predictions can be experimentally tested through interfacial impedance spectroscopy and depth-profiling analyses across different glass compositions.

4. Conclusion

In this work, we developed and applied an ML-FF to investigate the early-stage interfacial reactions and SEI formation between lithium metal and three representative glassy sulfide electrolytes: LiSiS, LiSiPS, and LiPS. Large-scale MD simulations reveal a unified picture of early-stage SEI evolution in these disordered systems, characterized by spontaneous formation of a Li₂S-rich amorphous interphase and growth that slows over time, consistent with transport-limited behavior over the simulated window. We show that electrolyte composition plays a central role in governing interfacial reactivity and growth rates. P-containing glasses exhibit faster SEI formation, consistent with more extensive P–S bond reduction, whereas Si-rich systems display slower kinetics and retain partially reduced Si–S units near the interface. These compositional effects lead to distinct interphase chemistries and transport properties. In LiPS, the development of sufficiently large Li₂S domains can, in some trajectories, trigger stochastic crystallization, yielding a P-containing Li₂S-rich phase with high Li⁺ mobility and a possible passivating capability. In contrast, Si-containing systems remain predominantly amorphous and form mixed ionic–electronic conducting interphases that are intrinsically less stable.

Our results further demonstrate that, despite their electronically insulating nature, both amorphous and partially crystalline SEIs remain permeable to Li⁺ transport, enabling continued electrochemical operation. The high ionic conductivity of the crystallized Li₂S-rich interphase arises from aliovalent P incorporation and vacancy formation, revealing the importance of defect chemistry in regulating interfacial performance. By systematically comparing glassy and crystalline sulfide interfaces, this study clarifies how structural disorder, network chemistry, and reaction kinetics collectively determine SEI morphology and stability. While thermodynamics drives all sulfide electrolytes toward similar reduction products, the resulting interphase structure and properties are governed by kinetic pathways and local atomic environments. These findings bridge the gap between idealized crystalline models and practical glassy electrolytes used in real devices.

Author contributions

Conceptualization & project administration: Q. A. Investigation and methodology: R. Z., K. L. Supervision: Q. A. Writing – original draft: R. Z. Writing – review & editing: all authors. Resources and funding acquisition: Q. A.

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: additional ML-FF validation, SEI thickness fitting, RDF and coordination-number analyses, Li|SEI interface evolution, and SEI thickness time-dependence model comparisons. See DOI: https://doi.org/10.1039/d6ta01253a.

Acknowledgements

This work was supported by the start-up grant from Iowa State University, and the simulations were performed at ISU High Performance Computing clusters.

References

1. J. Janek and W. G. Zeier, Nat. Energy, 2016, 1, 16141.
2. X. Zeng, M. Li, D. Abd El-Hady, W. Alshitari, A. S. Al-Bogami, J. Lu and K. Amine, Adv. Energy Mater., 2019, 9, 1900161.
3. N. Kamaya, K. Homma, Y. Yamakawa, M. Hirayama, R. Kanno, M. Yonemura, T. Kamiyama, Y. Kato, S. Hama, K. Kawamoto and A. Mitsui, Nat. Mater., 2011, 10, 682–686.
4. J. Lee, T. Lee, K. Char, K. J. Kim and J. W. Choi, Acc. Chem. Res., 2021, 54, 3390–3402.
5. N. Boaretto, I. Garbayo, S. Valiyaveettil-SobhanRaj, A. Quintela, C. Li, M. Casas-Cabanas and F. Aguesse, J. Power Sources, 2021, 502, 229919.
6. Y. Zhu, X. He and Y. Mo, ACS Appl. Mater. Interfaces, 2015, 7, 23685–23693.
7. Y. Xiao, Y. Wang, S.-H. Bo, J. C. Kim, L. J. Miara and G. Ceder, Nat. Rev. Mater., 2019, 5, 105–126.
8. N. L. Marana, S. Casassa, M. F. Sgroi, L. Maschio, F. Silveri, M. D'Amore and A. M. Ferrari, Langmuir, 2023, 39, 18797–18806.
9. M. D'Amore, M. Y. Yang, T. Das, A. M. Ferrari, M. M. Kim, R. Rocca, M. Sgroi, A. Fortunelli and W. A. Goddard, J. Phys. Chem. C, 2023, 127, 22880–22888.
10. Y. Zhang, W. Wang, J. Ju, S. Zhang, P. Xu, J. Lu, Y. Yang, T. Dong, L. Cui, S. Liu, L. Ni, Z. Lv, P. Han, J. Xu, S. Cui, F. Sun, Y. Jin and G. Cui, Adv. Energy Mater., 2026, 16, e06260.
11. Q. Wu, S. Xiong, F. Li and A. Matic, Batter. Supercaps, 2023, 6, e202300321.
12. H. Adenusi, G. A. Chass, S. Passerini, K. V. Tian and G. Chen, Adv. Energy Mater., 2023, 13, 2203307.
13. Z. Zhang, Y. Li, R. Xu, W. Zhou, Y. Li, S. T. Oyakhire, Y. Wu, J. Xu, H. Wang, Z. Yu, D. T. Boyle, W. Huang, Y. Ye, H. Chen, J. Wan, Z. Bao, W. Chiu and Y. Cui, Science, 2022, 375, 66–70.
14. F. Walther, R. Koerver, T. Fuchs, S. Ohno, J. Sann, M. Rohnke, W. G. Zeier and J. Janek, Chem. Mater., 2019, 31, 3745–3755.
15. R. Zhao, G. Hu, S. Kmiec, R. Gebhardt, A. Whale, J. Wheaton and S. W. Martin, ACS Appl. Mater. Interfaces, 2021, 13, 26841–26852.
16. V. Torres, P. Philipp, S. Kmiec and S. W. Martin, Chem. Mater., 2024, 36, 5521–5533.
17. A. M. Nolan, Y. Zhu, X. He, Q. Bai and Y. Mo, Joule, 2018, 2, 2016–2046.
18. L. E. Camacho-Forero and P. B. Balbuena, J. Power Sources, 2018, 396, 782–790.
19. M. Sun, W. Yuan, X. Zhang, Y. Tian and Z. Zhou, npj Comput. Mater., 2025, 11, 350.
20. O. T. Unke, S. Chmiela, H. E. Sauceda, M. Gastegger, I. Poltavsky, K. T. Schütt, A. Tkatchenko and K.-R. Müller, Chem. Rev., 2021, 121, 10142–10186.
21. F. Ren, Y. Wu, W. Zuo, W. Zhao, S. Pan, H. Lin, H. Yu, J. Lin, M. Lin, X. Yao, T. Brezesinski, Z. Gong and Y. Yang, Energy Environ. Sci., 2024, 17, 2743–2752.
22. G. Chaney, A. Golov, A. Van Roekeghem, J. Carrasco and N. Mingo, ACS Appl. Mater. Interfaces, 2024, 16, 24624–24630.
23. Kausthubharam, B. S. Vishnugopi, A. S. J. Alujjage, V. Premnath, W. S. Tang, J. A. Jeevarajan and P. P. Mukherjee, Chem. Rev., 2025, 126(1), 404–447.
24. R. Zhao, G. Hu, S. Kmiec, J. Wheaton, V. M. Torres and S. W. Martin, Batter. Supercaps, 2022, 5(11) DOI:10.1002/batt.202100356.
25. T. A. Yersak, C. Kang, J. R. Salvador, N. P. W. Pieczonka and M. Cai, Mater. Adv., 2022, 3, 3562–3570.
26. Z. Xu and Y. Xia, J. Mater. Chem. A, 2022, 10, 11854–11880.
27. H. Wang, L. Zhang, J. Han and W. E, Comput. Phys. Commun., 2018, 228, 178–184.
28. L. Zhang, J. Han, H. Wang, R. Car and W. E, Phys. Rev. Lett., 2018, 120, 143001.
29. J. Zeng, D. Zhang, A. Peng, X. Zhang, S. He, Y. Wang, X. Liu, H. Bi, Y. Li, C. Cai, C. Zhang, Y. Du, J.-X. Zhu, P. Mo, Z. Huang, Q. Zeng, S. Shi, X. Qin, Z. Yu, C. Luo, Y. Ding, Y.-P. Liu, R. Shi, Z. Wang, S. L. Bore, J. Chang, Z. Deng, Z. Ding, S. Han, W. Jiang, G. Ke, Z. Liu, D. Lu, K. Muraoka, H. Oliaei, A. K. Singh, H. Que, W. Xu, Z. Xu, Y.-B. Zhuang, J. Dai, T. J. Giese, W. Jia, B. Xu, D. M. York, L. Zhang and H. Wang, J. Chem. Theory Comput., 2025, 21, 4375–4385.
30. R. Zhou, K. Luo, S. W. Martin and Q. An, ACS Appl. Mater. Interfaces, 2024, 16, 18874–18887.
31. S. P. Ong, W. D. Richards, A. Jain, G. Hautier, M. Kocher, S. Cholia, D. Gunter, V. L. Chevrier, K. A. Persson and G. Ceder, Comput. Mater. Sci., 2013, 68, 314–319.
32. J. Zeng, X. Peng, Y.-B. Zhuang, H. Wang, F. Yuan, D. Zhang, R. Liu, Y. Wang, P. Tuo, Y. Zhang, Y. Chen, Y. Li, C. T. Nguyen, J. Huang, A. Peng, M. Rynik, W.-H. Xu, Z. Zhang, X.-Y. Zhou, T. Chen, J. Fan, W. Jiang, B. Li, D. Li, H. Li, W. Liang, R. Liao, L. Liu, C. Luo, L. Ward, K. Wan, J. Wang, P. Xiang, C. Zhang, J. Zhang, R. Zhou, J.-X. Zhu, L. Zhang and H. Wang, J. Chem. Inf. Model., 2025, 65, 11497–11504.
33. Y. Zhang, H. Wang, W. Chen, J. Zeng, L. Zhang, H. Wang and W. E, Comput. Phys. Commun., 2020, 253, 107206.
34. G. Kresse and J. Furthmüller, Comput. Mater. Sci., 1996, 6, 15–50.
35. G. Kresse and J. Furthmüller, Phys. Rev. B, 1996, 54, 11169–11186.
36. G. Kresse and J. Hafner, Phys. Rev. B, 1993, 47, 558–561.
37. J. P. Perdew, A. Ruzsinszky, G. I. Csonka, O. A. Vydrov, G. E. Scuseria, L. A. Constantin, X. Zhou and K. Burke, Phys. Rev. Lett., 2008, 100, 136406.
38. A. P. Thompson, H. M. Aktulga, R. Berger, D. S. Bolintineanu, W. M. Brown, P. S. Crozier, P. J. in ’t 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.
39. L. Martínez, R. Andrade, E. G. Birgin and J. M. Martínez, J. Comput. Chem., 2009, 30, 2157–2164.
40. A. Kato, H. Kowada, M. Deguchi, C. Hotehama, A. Hayashi and M. Tatsumisago, Solid State Ionics, 2018, 322, 1–4.
41. S. Wenzel, D. A. Weber, T. Leichtweiss, M. R. Busche, J. Sann and J. Janek, Solid State Ionics, 2016, 286, 24–33.
42. L. Sang, R. T. Haasch, A. A. Gewirth and R. G. Nuzzo, Chem. Mater., 2017, 29, 3029–3037.
43. C. Szczuka, B. Karasulu, M. F. Groh, F. N. Sayed, T. J. Sherman, J. D. Bocarsly, S. Vema, S. Menkin, S. P. Emge, A. J. Morris and C. P. Grey, J. Am. Chem. Soc., 2022, 144, 16350–16365.
44. Z. Zheng, K. Wang, Z. Zhu, Y. Dai, S. Wang, Z. Wu, J. Zhou and Y. Zhu, Energy Storage Mater., 2026, 84, 104804.
45. P. M. Attia, W. C. Chueh and S. J. Harris, J. Electrochem. Soc., 2020, 167, 090535.
46. C. D. Alt, N. U. C. B. Müller, L. M. Riegger, B. Aktekin, P. Minnmann, K. Peppler and J. Janek, Joule, 2024, 8, 2755–2776.

---