Atomic-level insights into the highly conductive lithium thio-phosphate solid electrolytes with exceptional stability against lithium metal

Abstract

Despite the wide range of emerging solid electrolytes with promising characteristics, such as high ionic conductivity, the inherent thermodynamic instability against lithium metal remains a significant challenge. We have previously introduced a new family of solid electrolytes based on thiophosphates with inherent stability against the metallic Li anode and high conductivity [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]. In this study, we employ density functional theory (DFT) together with ab initio molecular dynamics (AIMD) simulations to investigate the diffusion mechanisms and underlying factors contributing to the high ionic conductivity of these novel Li–P–S ternary electrolytes, including Li₇PS₂, Li₅PS, Li₈P₂S, and Li₁₁P₃S. Our findings reveal that these materials exhibit ionic conductivities comparable to the well-known superionic conductor, Li₇P₃S₁₁, positioning them as promising candidates for solid-state battery applications. Additionally, we show that, unlike Li₇P₃S₁₁ which forms a solid-electrolyte interphase, the novel Li–P–S ternaries exhibit remarkable stability against Li metal anode due to their unique Li₂S-like structural framework. The absence of a solid-electrolyte interphase layer is particularly significant, as it eliminates additional resistance at the electrolyte–anode interface, a common challenge in many solid-state battery systems. Our study not only highlights the suitability of these novel ternaries, particularly Li₅PS, as high-performance solid electrolytes but also underscores the importance of structural design in developing next-generation battery materials. The ability of these materials to maintain high ionic conductivity and stability over extended periods makes them ideal candidates for future solid-state lithium batteries, offering a pathway to safer, more efficient, and longer-lasting energy storage solutions.

---

1. Introduction

In recent years, the pursuit of safer and more efficient energy storage systems has prompted the exploration of alternatives to conventional lithium-ion batteries (LIBs).¹ Traditional LIBs, fuelled by liquid electrolytes, exhibit inherent challenges related to safety and energy density. Liquid electrolytes are typically composed of organic solvents, such as ethylene carbonate (EC) or dimethyl carbonate (DMC), which are highly flammable. In the event of a short circuit, overheating, or mechanical damage to the battery, these solvents can ignite, leading to fires or explosions.^2,3 Lithium metal has the highest theoretical capacity of any anode material, making it highly desirable for high-energy-density batteries.^4–6 However, to fully exploit this potential, the electrolyte must be stable in contact with lithium metal. Lithium metal is highly reactive, especially when in contact with conventional liquid electrolytes, which can lead to the formation of lithium dendrites.^7–11 These dendrites can grow through the electrolyte and eventually cause short circuits, leading to battery failure, overheating, and potential fires or explosions. The emergence of all-solid-state batteries (ASSBs), with solid-state electrolytes replacing the liquid counterparts, has marked a significant paradigm shift. A solid electrolyte (SE) may enable the use of lithium metal anodes, leading to batteries with significantly higher energy densities compared to those using other anode materials.^5,10,12–15 Therefore, the transition toward ASSBs represents a transformative leap forward because it not only potentially allows the use of pure or alloyed lithium anodes, enabling higher capacities, but also holds the promise of mitigating flammability concerns.

ASSBs exhibit superiority over conventional LIBs, particularly in high-temperature applications, offering reduced flammability risks and a diminished threat of thermal runaway.^16,17 Furthermore, the absence of bulk polarization effects, as Li⁺-ions exclusively conduct ionic charge, suggests the potential for faster (dis)charging in ASSBs compared to traditional LIBs.¹⁸ Despite these advantages, the commercialization of ASSBs faces challenges, primarily the absence of a highly conductive SE that is stable toward lithium metal and cathode materials, while accommodating mechanical stresses induced during battery cycling. Dendrite-free plating of lithium metal also remains an unresolved issue.

The exploration of promising SEs has witnessed a proliferation of interest in diverse materials. One class, represented by cation-substituted Li₇La₃Zr₂O₁₂ (LLZO)-derived compounds,¹⁹ offers high ionic conductivity and kinetic stability toward lithium metal but faces challenges related to mechanical stiffness and susceptibility to cracking as well as lithium dendrite formation particularly at high current densities.¹⁸ A second promising class involves materials containing sulphide ions, exhibiting extremely high ionic conductivity and stress accommodation.²⁰ In these materials, oxide anions are substituted by larger and more polarizable sulphide anions. This results in reduced Li-ion jump barriers which enhances their conductivity enormously.²¹ In fact, materials incorporating thiophosphate polyhedra as fundamental units exhibit remarkably high ionic conductivities comparable to those observed in liquid electrolytes. Studying glasses and crystalline phases within the pseudo-binary (Li₂S)_x–(P₂S₅)_1−x were one of the first attempts to utilize the presence of thiophosphate polyhedra.^22,23 Recent developments have highlighted that doped crystalline derivatives such as Li₁₀GeP₂S₁₂ (LGPS, σ = 12 mS cm⁻¹),²⁴ Li_6−xPS_5−xClBr_x (σ = 24 mS cm⁻¹),²⁵ and Li_9.54Si_1.74P_1.44S_11.7Cl_0.3 (σ = 25 mS cm⁻¹)²⁶ are some of the most conductive SEs to date. More recently, a third class, lithium-rich ternary phosphides emerged as potential SE candidates, featuring phosphide anions (P³⁻) with even greater polarizabilities than sulphide anions. These phases are structured around anionic TtP₄ tetrahedra (Tt = Al, Si, Ge, Ga, or Sn). Notably, Li₉AlP₄ demonstrates conductivities of up to 3 mS cm⁻¹ at room temperature.^27–33 Despite their appeal as SEs due to a high number of charge carriers and low density, these materials exhibit sensitivity to oxygen and moisture.

Although the newly introduced SEs offer several advantages, they share a significant drawback: they are inherently thermodynamically unstable when in contact with lithium metal, as demonstrated for metal/metalloid-containing oxides, sulphides, and thiophosphates.^34–37 For instance, thiophosphates degrade into a mixture of lithium sulphides and phosphides, driven by the reduction of phosphorus ions from a P⁵⁺ oxidation state to P³⁻ in Li₃P. The resulting lithium binary compounds (Li₂O, Li₂S, Li₃P, and LiX (X = a halide)) are suboptimal as electrolytes, creating a conductivity bottleneck at the SE//Li metal interface. Additionally, the reduction of other cations in the material may lead to short circuits via the formation of lithium intermetallics. For example, Ge⁴⁺ in LGPS³⁸ and phosphidogermanates can form electronically conductive germanides,^36,37 posing a risk of further decomposition and short circuits through the creation of mixed ionic-electronic conducting interphases.³⁸

Quantum chemical calculations employing density functional theory (DFT) have been widely used to investigate the structural, electronic, transport, and spectral properties of various Li-ion battery materials, particularly SEs.^39–41 Moreover, DFT is often combined with crystal structure prediction approaches, such as Ab Initio Random Structure Search (AIRSS), genetic algorithm (GA), USPEX, and CALYPSO, in the in silico discovery of novel functional materials.^42,43 While these tools have been applied to predict chemically-doped thiophosphate electrolytes, such as Li₃Y(PS₄)₂ (ref. 44) and Li_x(PS₄)_yX_z (X = Cl, Br, I),⁴⁵ the pure Li–P–S phase diagram has largely remained unexplored, with the exception of individual studies focusing exclusively on a few Li₂S–P₂S₅ phases.²²

Although novel SEs show promise, their instability with lithium metal underscores the need for a comprehensive approach focusing on both structural and thermodynamic properties, as well as synthesis methods to developing highly Li⁺-ion conductive SEs that are thermodynamically stable against Li metal. In our previous work,⁴⁶ we successfully explored the Li–P–S ternary phase space, employing a combination of high-throughput crystal structure predictions and solid-state synthesis (via ball milling) to isolate the most promising compositions, particularly those within the Li₃P–Li₂S tie line. We conducted a systematic characterization of the structural properties and Li-ion mobility of these materials using techniques such as X-ray and neutron diffraction, solid-state nuclear magnetic resonance spectroscopy (relaxometry), and electrochemical impedance spectroscopy. This led to the identification of several Li₃P–Li₂S metastable solid solutions, where the phases adopted a fluorite (Li₂S) structure with phosphorus substituting for sulphur and extra Li⁺-ions occupying octahedral voids. The combined analysis of experimental data and quantum-chemical calculations⁴⁶ reveals that four of the new ternary compounds—Li₇PS₂, Li₅PS, Li₈P₂S, and Li₁₁P₃S—not only exhibit high structural stability and ionic conductivity (hull distances of 76.3 meV, 38.5 meV, 12.4 meV, 45.4 meV, respectively) but also possess low activation barriers for Li⁺-ion transport (193 meV, 190 meV, 232 meV, 250 meV at room temperature, respectively). Moreover, each of these structures are shown to be insulators with a direct band gap of at least 1.4 eV making them perfect candidates as SEs if they can maintain thermodynamic stability at the interfaces with electrodes.⁴⁶ As shown in Fig. 1, these compounds adopt a Li₂S-like orthorhombic (slightly distorted cubic) structural framework, which suggests potential thermodynamic stability in the presence of lithium metal anodes similar to Li₂S.^47–49

Fig. 1 Supercells of highly conductive Li₃P–Li₂S solid solutions studied in this work; (a) Li₇P₂S, (b) Li₅PS, (c) Li₈P₂S and (d) Li₁₁P₃S.

In this work, we delve into the thermodynamic stability of these four Li–P–S ternaries (see Fig. 1) in contact with lithium metal anodes by means of DFT simulations. The results are compared with established Li–P–S electrolytes (Li₇P₃S₁₁ and γ-Li₃PS₄) and end members (Li₂S and Li₃P) of the tie-line (see Fig. S1† for their unit cells). The article is structured as follows: Section 2 details the methodology and computational approach employed, Section 3 presents the results, beginning with an examination of the bulk structures, followed by an analysis of the surface structures, and concluding with an evaluation of the interfaces and their properties, and finally, Section 4 offers a summary of our findings and outlines the implications for future research.

Our approach represents a significant advancement in the development of SEs by addressing critical challenges in the field, particularly the instability of SE materials against lithium metal anodes. The novel Li–P–S ternaries we study not only demonstrate exceptional stability against lithium anodes—comparable to the highly stable Li₂S—but also exhibit high ionic conductivity without the formation of a solid-electrolyte interphase (SEI), even at elevated temperatures. This absence of SEI formation eliminates internal resistance at the interface, maintaining the high conductivity observed in their bulk structures. Furthermore, our investigation into the underlying mechanisms reveals that phosphorus substitution for sulphur atoms stabilizes extra interstitial sites and coordinates additional lithium atoms, creating new conduction channels that further enhance ionic transport. This dual achievement of stability and conductivity presents a pivotal step forward, offering a robust solution to the longstanding issues in SE design and paving the way for more efficient and safer energy storage systems.

2. Computational methods

2.1. INTERFACER code

To investigate the stability, SEI layer formation and ionic conductivity at electrolyte–electrode interfaces, we built model interfaces between Li–P–S structures and lithium metal. To facilitate this, we have developed the INTERFACER code,⁵⁰ an automated interface generation tool designed to systematically explore various candidate interfaces with minimal strain, starting from predicted bulk structures. By searching through specified Miller indices, INTERFACER identifies energetically-relevant interfaces and, in conjunction with ab initio DFT codes (Vienna Ab Initio Simulation Package (VASP) or Cambridge serial total energy package (CASTEP)), determines the thermodynamically and mechanically stable configurations. The INTERFACER tool, coded in Python, also integrates parameter setup for ab initio simulations, streamlining the interface optimization process.

2.2. Parameters for relaxation and AIMD simulations

Plane-wave DFT electronic structure calculations were performed using VASP (v.6.3),^51–53 which is an implementation of periodic boundary conditions and the pseudopotential approximation. We employed the projector-augmented wave (PAW)^54,55 method jointly with the Perdew-Wang (PW91) version of the generalized gradient approximation (GGA) exchange–correlation potentials.⁵⁶ The atomic positions and lattice parameters were fully relaxed using conjugate-gradient method until all forces acting on atoms were smaller than 0.01 eV Å⁻¹. The Brillouin zone was sampled using a Γ-centred Monkhorst–Pack (MP) grid with a k-point spacing finer than 0.1 Å⁻¹. The normal accuracy setting (PREC = Normal) along with an electronic convergence criterion of 10⁻⁸ eV and a Gaussian smearing factor of 0.1 eV were adopted. A cut-off-energy for the planewave basis functions of 520 eV was adopted for the relaxation calculations, while the cut-off was lowered to 400 eV for the ab initio molecular dynamics (AIMD) simulations. To allow for unbiased ionic migration dynamics, no symmetry constraints were applied during the dynamic simulations. AIMD simulations were performed not only to analyse the integrity of the interface, i.e. the SEI formation, but also for estimating the Li-ion conductivity within the Li–P–S ternaries as well as through the interfaces under investigation. For AIMD simulations, we benefitted from the faster implementation in VASP designed for the Γ-only calculations, which proved useful for obtaining numerous adequately long trajectories to yield better analysis for diffusivity. All crystal structures were visualized with VESTA⁵⁷ and Ovito.⁵⁸

3. Results and discussion

3.1. Bulk conductivity: the effect of phosphorus

Efficient ion transport in SEs is essential for achieving high-performance ASSBs. Consequently, understanding the diffusion processes and optimising ionic conduction at the nanoscale is a critical focus of contemporary research. Many thiophosphate-based solid electrolytes have been systematically studied to investigate how their average structures, local bonding environments, and lattice dynamics govern ionic conductivity.^22,59–66 For example, Li₇P₃S₁₁-type sulphides, composed of both PS₄³⁻ and P₂S₇⁴⁻ units, exhibit high conductivities owing to their flexible framework and interconnected polyhedral network, which facilitate lithium diffusion pathways.^65,66 Recent experimental and computational studies have further highlighted that local structural heterogeneity and polyhedral disorder can critically influence transport properties beyond what is predicted by average crystallographic models.^22,64,65 Moreover, state-of-the-art neutron scattering and first-principles simulations have revealed that certain thiophosphates display liquid-like lithium dynamics in the solid state, with anharmonic vibrations and concerted ion-hopping mechanisms contributing to their superionic behavior.⁵⁹ These findings underscore the importance of incorporating both structural and dynamic considerations—at both average and local levels—when interpreting ion transport in thiophosphates and developing new solid electrolytes.

In our previous study,⁴⁶ we presented both computational and experimental data on the ionic conductivities and activation energies of newly synthesised and other commonly studied Li–P–S ternary systems. Our findings revealed that the new ternaries, particularly Li₇PS₂, Li₅PS, Li₈P₂S, and Li₁₁P₃S (see Fig. 1) exhibit exceptionally high conductivities and low activation energies, for which we are yet to clarify the underlying reasons, comparable to those of the superionic conductor Li₇P₃S₁₁. While traditional experimental techniques offer valuable insights into macroscopic ionic conductivity, they often lack the resolution necessary to elucidate the atomic-scale mechanisms that govern ion movement, particularly in complex solid-state systems. AIMD enables the detailed investigation of microscopic processes underlying ionic diffusivity, offering an atomistic perspective of ion transport that accounts for the intricate interactions within the material's crystal lattice and electronic structure in a fully quantum mechanical framework. These simulations are therefore indispensable for the design and optimisation of new SE materials with enhanced ionic conductivity, directly influencing the development of next-generation energy storage technologies. In this section, we provide an understanding of the fundamental mechanisms behind the superior diffusion properties of the four newly-developed ternaries through AIMD simulations.

The self-diffusivity for particle i, also referred to as the tracer diffusivity (D*), can be computed from the mean square displacement (MSD) of the particle during an AIMD simulation as:

where d is the dimensionality of the diffusion, which we use 3 in all cases presented here. The MSD of the N particles within a bulk structure is expressed as:

The bulk ionic conductivity of an electrolyte, such as a Li–P–S ternary, can be estimated using the Nernst–Einstein relation,⁶⁷ given by:

where n is the density of diffusing particle, e is the elementary charge, Z is the ionic charge, k is the Boltzmann constant, T is the temperature and D is the ion diffusivity, with the Haven ratio (H_R) assumed to be unity (indicating non-correlated ionic motion). Furthermore, the activation energy for a given electrolyte can be determined from an Arrhenius plot, which involves the exponential fit of conductivity versus inverse temperature over a range of temperatures.

Fig. 2a illustrates the MSD of lithium atoms for bulk-Li–P–S systems at 525 K. Li₂S and Li₃P are not included in this figure due to their almost negligible diffusivities at 525 K. The data reveals that while the lithium mobility in the superionic conductor Li₇P₃S₁₁ is higher than in all other systems, the mobilities in Li₅PS and Li₈P₂S are quite comparable. Notably, while Li₇P₃S₁₁ has the highest tracer diffusivity (D*), it is only about four times that of Li₅PS and Li₈P₂S, indicating that the difference in mobility is moderate. Additionally, all the novel ternary compounds demonstrate significantly higher lithium diffusivity compared to γ-Li₃PS₄ meanwhile Li₅PS and Li₈P₂S ternaries present higher lithium diffusivities than that of highly conductive β-Li₃PS₄ (see Table S1† for the list of diffusivities and the calculated error values^68,69 for the systems shown in Fig. 2a). These observations underscore the enhanced ionic transport properties of the newly developed ternary systems.

Fig. 2 Diffusion properties of bulk Li–P–S ternaries at 525 K. (a) MSD of Li–P–S ternaries. Atomistic structure (upper panel) and Li⁺-ion trajectory (lower panel) for (b) Li₂S and (c) Li₇PS₂. (d) Conductivity of Li–P–S systems as a function of Li₂S–Li₃P composition ratio. Tracer diffusivities for Li₇P₃S₁₁, Li₅PS and Li₈P₂S are shown in (a) in red, green and orange colours, respectively. The structures in (b) and (c) are depicted along [111] for ease of comparison to interface structures. The red, black and blue arrows indicate possible diffusion channels for Li⁺-ions along 8c–8c (T–T), 8c–4b–8c (T–O–T) and 4b–4b (O–O), respectively. T/O denotes the tetrahedral/octahedral sites.

From AIMD trajectories, one can also extract the Li-ion transport pathways to elucidate the underlying factors contributing to the conductivity. In the Li₂S structure, lithium atoms occupy the 8c tetrahedral sites, while sulphur atoms reside in the 4a octahedral sites, as illustrated in the upper panel of Fig. 2b. The 4b octahedral sites are energetically unfavourable for both sulphur and lithium atoms, rendering them unstable and vacant. Three primary ion migration paths can be anticipated in Li₂S and Li–P–S ternaries: (1) the 8c–8c (T–T) path as shown by red arrows, involving Li hopping between thermodynamically stable tetrahedral sites (i.e., vacancy hopping);^40,70 (2) the 8c–4b–8c (T–O–T) path as shown by black arrows, where Li jumps between regular and interstitial Li sites;^40,70 and (3) the 4b–4b (O–O) path as shown by blue arrows, which is only possible for Li–P–S ternaries as Li₂S lacks interstitials. In the lower panel of Fig. 2b, we depict the lithium ion trajectory during the whole simulation of 100 ps at 525 K with green dots, while the initial positions of sulphur atoms are given as reference. For pristine Li₂S, Li conduction does not occur at 525 K. Li mobility is rather limited and confined to the regular Li sites, primarily due to the full occupancy of the 8c sites in the model used. Although Li migration through the 8c–8c (T–T) path appears to be the most probable diffusion mechanism, introducing Li lattice vacancies at the 8c sites is not expected to significantly enhance Li mobility at low temperatures.^40,70 Moreover, T–O–T type migration is also absent, consistent with the reported high activation barrier (E_a = 0.39 eV).⁴⁰ However, at 1075 K, slightly above the previously reported superionic phase transition temperature (T ≈ 900 K),^70–72 Li ions begin to diffuse via the T–O–T path, which then becomes the predominant mode of conduction.

Next, we investigate the ionic conduction in the novel Li₂S–Li₃P solid solutions. These novel ternary compounds exhibit significantly higher ionic conductivities compared to Li₂S, primarily due to the incorporation of phosphorus atoms. In these ternary compounds, each phosphorus atom substitutes a sulphur atom, occupying a 4a site. Phosphorus atoms, with their 3⁻ oxidation state, exert stronger coulombic interactions with lithium ions compared to sulphur atoms with a 2⁻ oxidation state. This interaction not only stabilises the adjacent 4b sites, enabling the coordination of additional lithium atoms, as depicted in the upper panel of Fig. 2c, but also reduces the energy barrier for T–O migration. Consequently, diffusion channels primarily form around the phosphorus atoms and most lithium diffusion in these ternary compounds occurs in the vicinity of phosphorus atoms due to the increased availability of stable interstitial sites, as demonstrated for bulk-Li₇PS₂ structure shown in the lower panel of Fig. 2c. Substituting S atoms with P does not appear to alter the 3D Li migration pathways significantly, as both T–T (more visible in other ternaries) and T–O–T paths (but not O–O) are utilised in all ternaries. A clear distinction between the parent Li₂S and the Li₂S–Li₃P solutions is the significantly higher prevalence of T–O–T (8c–4b–8c) conduction channels even at low temperatures, making the T–O–T path the primary mode of conduction. This is likely facilitated by the additional interstitial Li ions. The diffusion process typically initiates via a knock-on mechanism, wherein a lithium ion occupying an interstitial 4b site facilitates the displacement of a lithium ion from a lattice 8c site into a neighbouring, higher-energy vacant 4b interstitial site, constituting the T–O segment of the T–O–T conduction pathway. Subsequently, the resulting vacancy at the 8c site is filled by another lithium ion migrating from a nearby 4b interstitial site, completing the O–T segment of the pathway. Due to this mechanism and the limited number of 4b sites the phosphorus-to-sulphur ratio becomes a critical factor in optimising ionic conduction. Fig. 2d shows the calculated lithium ion conductivity at 525 K as a function of phosphorus content. Li₂S and Li₃P, representing the extremes, exhibit negligible conductivity. The data indicate that conductivity peaks around a phosphorus content of 50% (Li₅PS) where the phosphorus atoms are anticipated to distribute throughout the electrolyte more homogenously. Insufficient phosphorus results in a low number of stabilised 4b sites, reducing the number of available hopping sites for lithium ions. Conversely, an excess of phosphorus leads to an over-coordination of lithium ions, decreasing the number of vacant, stable 4b sites and thus limiting ionic conduction.

Based on these findings, we propose that incorporating phosphorus sites into the Li₂S–Li₃P solutions enhances Li migration through two mechanisms: (1) by lowering the migration energy barrier between nearby 8c–4b (T–O) sites due to strong Coulomb interactions between Li⁺–P³⁻ ions, and (2) by coordinating additional Li ions at the higher-energy interstitial octahedral (4b) sites. The former mechanism aligns with previous reports on other superionic conductors.⁷³ Both factors facilitate migration along the 8c–4b–8c and 8c–8c paths, with the former being more prominent, while the 4b–4b path remains less populated.

Due to the fundamentally distinctive structural configurations of the superionic conductors Li₇P₃S₁₁ and Li₃PS₄ compared to Li₂S and novel ternary compounds, their lithium-ion conduction pathways are not directly comparable. In the case of Li₇P₃S₁₁, studies have shown that the dynamic flexibility of P₂S₇ ditetrahedral units facilitates Li⁺ diffusion by effectively reducing energy barriers and generating an open framework for ion transport. This mechanism enables lithium-ion mobility without relying on well-defined interstitial sites or vacancy-mediated diffusion.^{64–66,74,75} In contrast, the new Li–P–S solid solutions have a knock-on mechanism, which is facilitated by the additional Li atoms located at the interstitials that are coordinated by the phosphorous atoms.

3.2. Surface formation

Surface formation energy is a crucial parameter in materials science, quantifying the stability of a surface relative to the bulk material. The derivation of surface formation energy requires modelling of the bulk and surface (slab) structures and calculating their total energies separately.^76,77 By comparing these energies, the surface formation energy (γ_s) can be derived as the energy difference between the surface and the corresponding bulk per unit area as follows:Here, E_slab is the total energy of the slab, E_bulk is the energy of the bulk crystal per formula unit, n is the number of formula units in the slab model and 2A is the surface area of the two surfaces of the slab. The accuracy of this formula relies on the slab being symmetric with respect to its central plane—meaning the atomic structure, composition, and arrangement below the top surface must mirror those above the bottom surface—while also preserving the stoichiometric ratio of the constituent elements. However, in instances where keeping this mirror symmetry necessitates the breach of stoichiometry or where imposing the symmetry proves unattainable, alternative considerations arise. In scenarios where breaking stoichiometry is necessary, the calculation of surface formation energy can be achieved by adapting the formula to incorporate chemical potentials pertaining to individual elements or a predetermined set of molecules. This modification is essential to accommodate the energy contributions introduced by additional/missing atoms, and it can be expressed as follows:where j is the number of different kinds of the extra elements, c_i and µ_i are the number and chemical potential of the extra element of kind i, respectively. The inherent issue with this methodology lies in the non-uniqueness of the chemical potential values assigned to specific elements. Conventionally, the determination of chemical potentials involves evaluating the energies associated with a range of constituent subsystems that collectively form the overarching structure. The selection of these subsystems can vary, and there is no consensus on the optimal approach. For instance, Xie et al.⁷⁸ adopt a phase diagram-based approach, selecting subsystems from the end members of the tie line in which the final product lies, while Canepa et al.⁷⁹ and Gao et al.⁸⁰ utilize the decomposition products as subsystems for chemical potential calculations (decomposition approach). Another strategy involves using the bulk forms of the elements (elemental approach), although for some cases this may yield high inaccuracy due to the presence of numerous intermediate structures between bulk elements and the final product. Moreover, the chemical potentials for gaseous species can vary significantly depending on the pressure and temperature, which results in the surface formation energy spanning a wide range of values.^81,82 Alternatively, starting materials for the final structure can serve as the basis for subsystem selection (starting material approach). Karasulu et al.⁸³ demonstrated in their investigation of the LLZO system that the last two approaches indeed produce different values for chemical potentials and formation energies, albeit without altering the ordering of formation energies.

In our study, we assess the chemical potentials associated with each element (Li, S and P) using the decomposition approach, analysing energies inherent to bulk lithium, lithium sulphide (Li₂S), and lithium phosphide (Li₃P), or alternatively, employing the starting material approach, analysing lithium sulphide, lithium phosphide, and phosphorus pentasulfide (P₂S₅), as well as the elemental approach, analysing bulk lithium (bcc-Li), bulk sulphur, and bulk phosphorus (black phosphorus). Each distinct configuration of these subsystems yields disparate chemical potentials, necessitating comparison with experimental data to validate their accuracy. Due to the lack of empirical methods to check the validity of these diverse set of chemical potentials which leads to a substantial ambiguity and the significant disparities observed in the formation energies of the structures studied in this work when employing chemical potentials derived from the aforementioned set of subsystems (see ESI, Tables S2 and S3†), we have opted to forgo the incorporation of chemical potential considerations and instead maintain consistent stoichiometry across all scenarios (even if it results in symmetry-breaking in some cases) similar to the approach of Haruyama et al.⁸⁴

As previously stated, the determination of surface formation energy necessitates the utilization of symmetric surfaces, a requirement that may not be generally achievable across all Miller indices. Indeed, the presence of symmetric surfaces within the systems under investigation is notably scarce. Consequently, we have devised a methodological approach guided by a predefined set of principles governing the generation of slabs, with particular emphasis on prioritizing specific surface characteristics.

The set of principles can be outlined as follows:

I. Ensure the stoichiometric integrity: This principle arises from the inherent non-uniqueness of chemical potentials and the lack of testing mechanisms, where outcomes are intricately tied to the chosen parameters.

II. Where feasible, create symmetric surfaces: Surface symmetry is imperative for the application of eqn (4), facilitating the precise determination of surface formation energies. The slabs featuring symmetric surfaces are denoted by the symbol “s” within Table 1.

Table 1 Calculated surface formation energies of the studied material systems using various Miller indices. Symmetry of the slab surfaces: s = symmetric, q = quasi-symmetric, n = fully non-symmetric

γ s (J m^−2)	Li	Li3P	γ-Li3PS4	Li7P3S11	Li2S	Li7PS2	Li5PS	Li8P2S	Li11P3S
(001)	0.46 ^s	0.48 ^s	0.45^q	0.25^s	0.88^s	0.65^q	0.65^s	0.65^q	0.71^q
(010)			0.21^s	0.18^s		0.67^q	0.67^s	0.66^q	0.68^q
(100)		0.51^s	0.34^s	0.09 ^s		0.66^q	0.61^s	0.65^q	0.74^q
(011)	0.50^s	0.67^s	0.41^q	0.12^s	0.50^s	0.51^q	0.56^s	0.59^q	0.65^q
(011̄)			0.40^n	0.17^s		0.50^q	0.56^s	0.62^q	0.64^q
(101)			0.27^q	0.12^s		0.51^q	0.59^s	0.57^q	0.69^q
(101̄)			0.46^n	0.15^s		0.55^q	0.59^s	0.58^q	0.59^q
(110)		0.73^s	0.18 ^s	0.20^s		0.53^q	0.62^s	0.60^q	0.66^q
(11̄0)			0.18^s	0.20^s		0.55^q	0.62^s	0.57^q	0.65^q
(111)	0.54^s	0.85^s	0.33^n	0.18^s	0.33 ^s	0.34 ^q	0.43 ^s	0.46^q	0.56^q
(111̄)			0.29^n	0.16^s		0.36^q	0.43^s	0.42 ^q	0.54 ^q
(11̄1)			0.36^n	0.17^s		0.34^q	0.43^s	0.43^q	0.58^q
(1̄11)			0.50^n	0.13^s		0.37^q	0.43^s	0.44^q	0.56^q
(112̄)				0.22^s
(102̄)				0.21^s
(210)	0.51^s		0.54^s
(211)	0.54^s		0.29^n

---

III. In case symmetric surfaces are unattainable, form quasi symmetric surfaces: We focus on aligning the surface atoms while disregarding those situated beneath. This approach aims to maximize symmetry within the outermost layers (while minimizing the net dipole moment), thereby enabling the utilization of eqn (4) under the assumption of comparable formation energies between the two surfaces. The slabs featuring quasi-symmetric surfaces are denoted by the symbol “q” within Table 1.

IV. Maintain the integrity of the S–P bonds: While not a strict rule, preserving the integrity of the S–P bonds in PS₄ tetrahedra or P₂S₇ (if any) during surface formation prevents substantial energy penalties.

V. In cases where neither symmetric nor quasi-symmetric surfaces can be achieved, the surfaces are constructed while adhering solely to stoichiometric considerations. For the structures in this category, the results will only be important if the calculated formation energy for the fully non-symmetric slab is found to be the lowest among all Miller indices indicating the need for subsequent investigation into the corresponding structure. Fully non-symmetric slabs are denoted as “n” within Table 1.

We systematically studied the free surface formation for bcc-Li, Li₃P, γ-Li₃PS₄, Li₇P₃S₁₁, Li₂S and the new ternaries (Li₇PS₂, Li₅PS, Li₈P₂S, Li₁₁P₃S) for all Miller indices up to 1 (i.e., −1, 0, and 1). For some cases, we further analysed specific Miller indices to compare with the literature. The results are presented in Table 1 where we show the lowest formation energies in bold font and the corresponding surfaces are depicted in Fig. S2.† (001) surface in bcc-Li (symmetrically equivalent to (100) surface) and (001) surface in Li₃P present the lowest formation energies of 0.46 J m⁻² and 0.48 J m⁻², respectively, which are in perfect agreement with previous studies.^85–87 Our results indicate that for γ-Li₃PS₄, (110) surface has the lowest surface energy, at 0.18 J m⁻². The surfaces of γ-Li₃PS₄ have not been extensively explored in the literature. Lepley et al.⁸⁸ reported a surface formation energy of 0.32 J m⁻² for the (010) surface though our calculations yield a lower formation energy of 0.21 J m⁻² for the same Miller index. Discrepancies in energy assessments may stem from the myriad of conceivable surfaces for Li₃PS₄, rendering comprehensive analysis impractical. For instance, the (010) surface presents at least three distinct cuts, the fewest among all Miller indices, further compounded by the lack of symmetry in many slabs, introducing additional ambiguity to model selection and ensuing energy calculations. For Li₇P₃S₁₁, (100) surface yields the lowest surface formation energy of 0.09 J m⁻², consistent with prior literature.⁸⁹ Due to the cubic symmetry of Li₂S, we focus on the formation of surfaces corresponding to three Miller indices: (100), (110), and (111). Although the cubic symmetry is disrupted in the new ternary compounds, their structural framework remains analogous to Li₂S. Consequently, their surfaces can be categorized into three distinct families of Miller indices (100), (110), and (111), as illustrated in Fig. S3.† These families can easily be identified by their structural motifs when viewed from above the surface. For example, in (100) family, the anions form a square shape (see the black dashed lines in Fig. S3† top view), whereas, in (110) and (111) families, a parallelogram and a triangular shape emerges, respectively. Specifically, the (100) family includes (010) and (001) surfaces; the (110) family includes (11̄0), (101), (101̄), (011), and (011̄) surfaces; and the (111) family includes (111̄), (11̄1) and (1̄11) surfaces. Although the surfaces for different Miller indices within each family are not equivalent, they still exhibit similar surface formation energies (Table 1). The (111) surface of Li₂S demonstrates the lowest formation energy at 0.33 J m⁻², which aligns well with values reported in the literature.^78,85,87 The red circles and lines on the far left of Fig. S3† represent the Li atoms and Li–S bonds lost during the cleavage to create the surface. Notably, for the (111) surface, sulphur atoms lose only a single Li bond, whereas for the (100) and (110) surfaces, two Li bonds are broken. We note that the breaking of fewer bonds during surface formation leads to a smaller increase in total energy, which in turn results in lower formation energy. Given the similar framework to Li₂S, the new ternary compounds also exhibit the lowest surface formation energies for (111) type surfaces. Specifically, for Li₇PS₂ and Li₅PS, (111) surfaces have the lowest formation energies of 0.34 J m⁻² and 0.43 J m⁻², respectively while, the lowest formation energy for Li₈P₂S and Li₁₁P₃S is observed for (111̄) surfaces, with values of 0.42 J m⁻² and 0.54 J m⁻², respectively. Notably, the surface formation energy increases with the phosphorus content. In these systems, each phosphorus ion can coordinate additional lithium ions. Consequently, when the bulk structure is cleaved to form a surface, the phosphorus ions at the surface lose more lithium coordination than sulphur ions. This implies that with higher phosphorus content, more Li–P bonds are broken upon surface formation, thereby increasing the surface formation energy.

The Wulff shape, is a geometric method utilized to depict the equilibrium morphology of a crystal.^90,91 In Fig. 3, we present the Wulff shapes derived from the surface formation energies listed in Table 1. The shape for bcc-Li aligns perfectly with the literature,⁸⁶ showing (001) and (110) surfaces contributing 39% and 34% of the total surface area, respectively. Contributions from (111), (210), and (211) are approximately 9% each. The shape for Li₃P shows some differences from what Canepa et al.⁸⁵ obtained as they also investigate Miller indices higher than 1. The Wulff shape for Li₃P is composed of (100) and (001) type surfaces, with 65% and 35% surface areas, respectively. (110) and (101) surfaces dominate the Wulff shape of γ-Li₃PS₄ by contributing 58% and 33% to the surface area, respectively. For Li₇P₃S₁₁, the surface area contributions are distributed among many surfaces, with (100) surface contributing 31%, (011) surface 19%, (101) surface 16%, and (111̄) surface 13%. For Li₂S, the Wulff shape exclusively contains the (111) surface. In new ternaries, similar to Li₂S, the Wulff shape is dominated by the (111) family surfaces. However, as the phosphorus content increases, other surface families begin to appear, resulting in a less anisotropic shape. This is attributed to the rapid increase in surface formation energy for the (111) family with increasing phosphorus content. For Li₁₁P₃S, the Wulff shape comprises (111), (100), and (110) surfaces, contributing 78%, 14%, and 8% to the total area, respectively.

Fig. 3 Wulff Shapes for all of the studied systems.

3.3. Interface stability at 0 K

Upon elucidating the surface formation energies, we employed the INTERFACER⁵⁰ computational framework to construct interfaces between Li(100) and Li–P–S surfaces (the workflow of the methodology is briefly explained in ESI†) ensuring minimal surface formation energies (values in bold font in Table 1). Interface formation energy for the interface between materials a and b can be calculated as;where E_ab is the total energy of the interface system, n_a (n_b) denotes the number of formula units of material a (b), with corresponding bulk energy per formula unit E_a(E_b), and A represents the surface area of the interface. The division by 2A accounts for the presence of two interfaces in the sandwich model (no vacuum padding between surface slabs).

Due to the imposed boundary conditions in the simulation system, a lattice strain is artificially introduced while aligning the lattices of materials a and b. However, in real materials, such boundary conditions are absent, potentially resulting in negligible strain. Consequently, the calculated interface formation energy is usually overestimated, encapsulating an energy component attributed to in-plane strain. The actual interface energy (γ_ab^act) therefore exists within upper and lower bounds represented as γ_ab^min < γ_ab^act < γ_ab^max, where γ_ab^max = and γ_ab^min denotes the minimal interface energy devoid of strain energy (E^str), calculated as;

γ_ab^min = γ_ab − E^str (7)

The strain energy can be determined through various methods^47,48,92,93 but we introduce only two of them here. The first method involves subtracting the energy of relaxed slabs from that of the strained slabs, which we refer to as “Comparison” method. Adopted from Ferrari et al.,⁴⁷ we allow interface structure to relax to its minimum energy configuration, inducing strain in both slabs. Thus, the strain energy for each slab can be calculated as follows;

where E_slab-a^strained is the energy of “slab a” under strain and E_slab-a^relaxed is the energy of the slab when the in-plane lattice vectors are allowed to relax. Since this method uses slabs, it ignores any contribution that may arise due to the presence of the interface. The second method, proposed by Lepley et al.,⁴⁸ incorporates the contribution of the interface on strain, which we call “Extrapolation” method. It calculates the in-plane strain at the interface between materials a and b as:

γ_ab = γ_ab^lim + n_bc (9)

where γ_ab^lim is the interface energy in the coherent limit, n_b is the number of formula units of material b and c is a constant that is dependent on the material and the applied stress. Eqn (9) assumes that the lattices of the interface systems are fixed to the bulk values of material a and only material b is under strain. However, since our interface simulations involve relaxation of the entire system rather than fixing the lattices to bulk values of one material, eqn (9) is modified to include strain in both materials as:

γ_ab = γ_ab^lim + n_ac_a + n_bc_b (10)

Here, plotting γ_ab against n_a (n_b) yields a straight line with slope c_a (c_b) and intercept γ_ab^lim + n_bc_b (γ_ab^lim + n_ac_a) (see Fig. S4a†). Consequently, the strain values for each material E_a^str = n_ac_a (E_b^str = n_bc_b) is calculated separately, and the total strain (E_a^str) is determined as the sum of individual strains. γ_ab^lim can then be calculated by inserting individual strain values into eqn (10). The difference between the strain energies from the two methodologies γ^lnt = γ_ab^lim − γ^min describes the contribution of the interface to strain energy.

A significant limitation of the extrapolation method is the rapid escalation in the number of atoms as n is increased for the Li–P–S systems. For the majority of the studied interfaces, we employ the minimal n values for Li–P–S slabs, which are already considerably thick. Consequently, this renders the extrapolation approach impractical for our purposes. Nevertheless, we conducted a series of tests using both approaches to evaluate the differential impacts between the comparison and the extrapolation (see Fig. S4b†) methods for two exemplary systems: Li₂S(111)//Li(100) and Li₅PS(111)//Li(100). The results are summarized in Table 2. The strain values obtained from both methodologies are of the same order, and the effect of the interface (γ^int) is less than 10%. Hence, we applied the comparison method to the remainder of the interfaces studied.

Table 2 Calculated strain energies in J m⁻² with comparison (eqn (8)) and extrapolation (eqn (10)) methods

(J m^−2)	Li2S(111)//Li(100)		Li5PS(111)//Li(100)
	Comparison	Extrapolation	Comparison	Extrapolation
E Li ^str	2.9 × 10^−2	3.2 × 10^−2	6.1 × 10^−2	5.6 × 10^−2
E electrolyte ^str	5.4 × 10^−2	4.6 × 10^−2	8.8 × 10^−2	8.0 × 10^−2
E ^str	8.3 × 10^−2	7.9 × 10^−2	0.15	0.14
γ ^min	0.31	—	0.14	—
γ ^lim	—	0.32	—	0.15
γ ^lnt	0.01 (1.5%)		0.01 (7.5%)
c Li	—	1.3 × 10^−4	—	3.8 × 10^−4
c electrolyte	—	6.8 × 10^−4	—	1.7 × 10^−3

---

The work of adhesion (W_adh) represents the energy required to separate two adherent surfaces to an infinite distance and quantifies the mechanical stability of an interface. In practical terms, W_adh is defined as the energy necessary to separate the two slabs and can be formulated as follows,

We define it as the minimum work of adhesion since it inherently contains the strain energy which lowers the value. To get the maximum work of adhesion, the strain energy should be added to this value (W_adh^max = W_adh^min + E^str). Similar to the interface formation energy, the actual work of adhesion lies between W_adh^min and W_adh^max. Typically, low (positive value) interfacial energy and high (positive value) work of adhesion are indicative of strong interfacial stability.

The calculated values for the interface formation energy, work of adhesion and strain energy are summarized in Table 3. The calculated maximum and minimum formation energies for Li₃P(001)//Li(100) and Li₂S(111)//Li(100) interfaces are γ_Li3P//Li^max = 0.37 J m⁻², γ_Li3P//Li^min = 0.35 J m⁻², γ_Li2S//Li^max = 0.40 J m⁻² and γ_Li2S//Li^min = 0.31 J m⁻², respectively, which are in agreement with the literature.⁸⁵ The interface formation energies for new ternaries are even lower than Li₂S and Li₃P, indicating higher stability. Interface formation energies are typically positive because the formation of interfaces involves creating surfaces from the bulk forms of two systems, which are inherently lower in energy. However, the interfaces of Li₃PS₄ and Li₇P₃S₁₁ with Li metal exhibit negative interface formation energies. This anomaly indicates the presence of chemical activities at the interface that lowers the total energy of the interface system. We present the Li₅PS//Li and Li₇P₃S₁₁//Li interfaces as examples of inactive and active interfaces, respectively, in Fig. 4 (See Fig. S5† for all other structures). The dashed lines indicate the position of the interface before relaxation. Notably, the PS₄ tetrahedra near the interface breaks up and sulphur atoms diffuse across the Li₇P₃S₁₁//Li interface to coordinate with lithium atoms and form Li₂S-like structures as indicated by the oval (Fig. 4b, bottom panel), whereas no diffusion of sulphur or phosphorus atoms is observed across the Li₅PS//Li interface. The S–P bond breaking and sulphur atom diffusion results in the undercoordinated phosphorus atoms bond with lithium atoms in the vicinity. These changes in the bond formation from S–P to Li–P and Li–S significantly reduce the total energy of the system, leading to negative formation energy. Other works that study the apparent reactivity of Li₃PS₄//Li and Li₇P₃S₁₁//Li interfaces agree on the formation of a similar Li₂S-like buffer layer.^{47,48,88,94,95}

Table 3 Calculated interface formation energies (γ_ab^max, γ_ab^min), adhesive energies (W_adh^max, W_adh^min) and the strain energy (from the comparison method) for Li–P–S systems against Li(100), all energies in J m⁻²

Interface	γ ab ^max	γ ab ^min	W adh ^min	W adh ^max	E ^str
Li3P(001)//Li(100)	0.37	0.35	0.56	0.58	2.0 × 10^−2
γ-Li3PS4(110)//Li(100)	−0.90	−1.11	1.41	1.63	0.22
Li7P3S11(100)//Li(100)	−1.87	−2.04	2.31	2.48	0.17
Li2S(111)//Li(100)	0.40	0.31	0.42	0.50	8.3 × 10^−2
Li7PS2(111)//Li(100)	0.27	0.25	0.50	0.52	1.9 × 10^−2
Li5PS(111)//Li(100)	0.28	0.14	0.49	0.64	0.15
Li8P2S(111)//Li(100)	0.25	0.22	0.56	0.60	3.2 × 10^−2
Li11P3S(111)//Li(100)	0.28	0.15	0.58	0.74	0.13

---

Fig. 4 The initial (upper) and the final (lower) states of (a) chemically inactive Li₅PS//Li, and (b) chemically active Li₇P₃S₁₁//Li interface structures. The dashed black line indicates the initial position of the interface. The black oval highlights the sulphur diffusion into lithium metal.

The calculated maximum and minimum work of adhesion for Li₃P(001)//Li(100) and Li₂S(111)//Li(100) interfaces are W_Li3P//Li^max = 0.58 J m⁻², W_Li3P//Li^min = 0.56 J m⁻², W_Li2S//Li^max = 0.50 J m⁻² and W_Li2S//Li^min = 0.42 J m⁻², respectively, which are consistent with the values reported in the literature.^85,87 The work of adhesion for the new ternary compounds with lithium is higher than that for Li₂S and Li₃P, indicating greater stability. Furthermore, the phosphorus content has an increasing effect on the work of adhesion, similar to the trend observed for surface formation energy. Having more phosphorus atoms, results in the formation of more bonds at the interface, making it more difficult to separate the surfaces. In the cases of Li₃PS₄//Li and Li₇P₃S₁₁//Li interfaces, the chemical reactions at the interface lead to a significant reduction in energy, yielding much higher work of adhesion compared to the other interfaces.

3.4. Dynamic interface stability at finite temperature

In our study, we employed AIMD simulations at 525 K and 1075 K to investigate the chemical interactions and stability of interfaces between Li-metal anode and the electrolytes as well as the effects of the interface on diffusion properties of Li ions at finite temperatures, with a particular emphasis on the formation of SEI layers. AIMD simulations provide high accuracy by solving quantum mechanical equations without relying on predefined potentials. Despite the computational expense limiting AIMD calculations to timescales of hundreds of picoseconds and system sizes to approximately 1000 atoms or fewer, these simulations furnish sufficient data to facilitate comparative analyses of different interfaces an provide critical insights into the structural integrity, dynamical properties, and the mechanisms governing SEI layer formation at the interface. In this section, we focus on the chemical interactions and the dynamical stability of the interfaces.

3.4.1. Stability of established Li–P–S electrolytes. First, we systematically investigated the time evolution of the intricate interfaces formed between lithium anode and commonly-studied lithium thiophosphate electrolytes, namely Li₂S, Li₃P, γ-Li₃PS₄, and Li₇P₃S₁₁. Since we found the behaviour of Li₃P//Li and γ-Li₃PS₄//Li in terms of dynamical stability and conductivity is very similar to that of Li₂S//Li and Li₇P₃S₁₁//Li, respectively, we will only focus on the latter two cases (the results for the former two cases can be found in Fig. S5–S8†). These simulations serve as a benchmark for understanding the behaviour of Li anode interfaces against the novel ternary systems. We present the evolution of the systems at various time steps for 525 K in Fig. 5a and b. The dashed black lines indicate the initial positions of the interfaces (at t = 0). Both Li₂S//Li and Li₃P//Li interfaces preserve their original structures over time even at an elevated temperature of 1075 K. Conversely, at the γ-Li₃PS₄//Li and Li₇P₃S₁₁//Li interfaces, diffusion of sulphur and phosphorus atoms into the lithium metal is observed, indicating instability akin to that at 0 K (i.e. geometry optimisations, see Fig. 4b). This instability intensifies and the degradation accelerates at higher temperatures.

Fig. 5 Evolution of the microstructure for (a) Li₂S//Li and (b) Li₇P₃S₁₁//Li interface structures at 525 K. The positional histogram of P and S atoms perpendicular to the interface for (c) Li₂S//Li and (d) Li₇P₃S₁₁//Li interface structures. The pair distributions of certain atomic pairs for (e) Li₂S//Li and (f) Li₇P₃S₁₁//Li interface structures. Black dashed lines in (a) and (b) indicate the initial positions of the interface.

The positional histograms (perpendicular to the interface – along the c-axis) for phosphorus and sulphur atoms, presented in Fig. 5c and d, elucidate the instability of the Li₇P₃S₁₁//Li interfaces, while simultaneously demonstrating the structural integrity of the Li₂S framework. In this stable framework, the S atoms wobble around their initial positions, retaining their original peak position, shape, and height in the histogram. In contrast, the P and S atoms in the highly-conductive Li₇P₃S₁₁ electrolyte become widely dispersed.

A similar conclusion can be drawn from the pair distribution functions (PDFs) for various atom pairs in these systems, as shown in Fig. 5e and f. In Li₂S//Li system, the separation distances between Li–S and S–S pairs remain relatively unchanged over time. However, in Li₇P₃S₁₁//Li system, the initial peak diminishes due to the dissolution of PS₄ tetrahedra at the interface, which decreases the number of P–S bonds.

Fig. 6a–c illustrate the distribution of atomic coordination numbers for various atoms at different times, using a coordination radius of 3 Å. While there is no significant change in Li₂S//Li system, in Li₇P₃S₁₁//Li system, the coordination number with S atoms remains relatively stable, whereas it increases for P atoms. This is due to the aforementioned dissolution of PS₄ units to form new Li–P bonds which increases the coordination number of P atoms from around 7 to 9.

Fig. 6 The coordination numbers for (a) S atoms in Li₂S//Li, (b) P and (c) S atoms in Li₇P₃S₁₁//Li interface structures. (d) The microstructure of Li₇P₃S₁₁//Li system after 20 ps at 525 K. The microstructures around P and S atoms within two interface subregions (decomposing and crystalline regions) as well as bulk Li₂S and bulk Li₃P are displayed for comparison. (e) The SEI layer thickness for γ-Li₃PS₄//Li and Li₇P₃S₁₁//Li as a function of time at 525 K. A significant drop in the SEI growth rate is observed at 83% and 85% consumption of P and S atoms in the Li₃PS₄//Li and Li₇P₃S₁₁//Li interface systems, respectively.

To gain a deeper understanding of the microstructural changes, we further analysed the Li₇P₃S₁₁//Li system. Given the interface instability, already at t = 20 ps, the interface region could be divided into two subregions, as shown in Fig. 6d: the decomposing region, also referred to as the SEI layer, and the crystalline region, which resembles the Li₇P₃S₁₁ crystal structure. Detailed examination of randomly-selected S and P atoms from these regions reveals significant differences in their coordination. Notably, P atoms, which are initially surrounded by S atoms in the crystalline region, become solely coordinated by Li atoms in the decomposing region. The coordination numbers of Li for P and S atoms in the decomposing region correspond well with those in the bulk Li₂S and bulk Li₃P. Therefore, a Li₂S-like layer where P atoms occasionally substitute S atoms forms, as suggested by other studies.^{47,48,88,94,95} With a much longer simulation on a much larger system, it might even be possible to observe the formation of separate Li₂S and Li₃P regions. This is, however, limited by the high cost of AIMD simulations, calling for other approaches, like machine-learned interatomic potentials. The decomposition reaction for the γ-Li₃PS₄ and Li₇P₃S₁₁ molecules are therefore as follows:

Li₃PS₄ + 8Li → 4Li₂S + Li₃P + 11.4 eV

and

Li₇P₃S₁₁ + 24Li → 11Li₂S + 3Li₃P + 35 eV

Here the reaction energies are calculated using the formation energies of the relevant species. These reactions result in a significant energy release due to the reduction of P atoms from 5⁺ to 3⁻.

By analysing the microstructure at each time step, we derived an estimation of the SEI layer thickness as a function of time, as depicted in Fig. 6e. Our model assumes that phosphorus and sulphur atoms coordinated solely with lithium atoms form Li₃P- and Li₂S-like structures, respectively. The SEI layer thickness was then calculated based on the relaxed bulk volumes of Li₃P and Li₂S as follows:

where d_SEI is the SEI layer thickness, n_P and n_S are the number of P and S atoms that are solely coordinated with Li atoms (coordination cut-off distance is 3 Å), respectively, V_Li3P and V_Li2S are the volumes per formula unit for relaxed Li₃P and Li₂S unit cells and A is the area of the interface. The SEI layer thickness was then calculated based on the relaxed bulk volumes of Li₃P and Li₂S. The analysis indicates that the initial structures, which are relaxed configurations at 0 K, already exhibit the presence of an SEI layer. During the time evolution, SEI layer formation follows a square root dependence on time for both the γ-Li₃PS₄//Li and Li₇P₃S₁₁//Li systems. After approximately 83-85% of the P and S atoms are consumed, the SEI growth rate decreases in both systems. The observed square root dependence of SEI layer thickness on time may be attributed to the already formed SEI layer hindering further growth, while the subsequent deceleration is likely a consequence of simulation limitations, specifically the depletion of available P and S atoms necessary for continued SEI layer formation at the initial rate. The higher rate of SEI formation in the Li₇P₃S₁₁//Li system compared to the γ-Li₃PS₄//Li system can be attributed to several factors. First, γ-Li₃PS₄ is structurally more stable than Li₇P₃S₁₁ in its bulk form,^22,89,96 which may reduce the propensity for SEI formation in the former. Second, the microstructures at the interfaces differ between the two systems. It is important to note that even within the same material system the SEI formation rates can vary depending on the Miller indices of the interfaces. Third, the systems are subjected to different levels of strain, which could influence the SEI formation process. Finally, the mechanisms of SEI formation differ between the two systems: in the γ-Li₃PS₄//Li system, SEI formation is primarily driven by the dissolution of PS₄³⁻ tetrahedral units, whereas in the Li₇P₃S₁₁//Li system, both PS₄³⁻ and P₂S₇⁴⁻ units contribute to the SEI layer formation.

3.4.2. Stability of novel Li–P–S compounds. In this section, we examine the stability and SEI layer formation at the interfaces between novel ternary compounds and lithium metal using AIMD simulations. Given that all the novel ternary-Li interfaces exhibit similar behaviour, we focus on the system with the most conductive ternary, namely Li₅PS//Li system, as a representative case and compare the results with the benchmark interfaces, viz. Li₂S//Li and Li₇P₃S₁₁//Li (The results for other ternary compounds can be found in Fig. S5–S8†).

Fig. 7a illustrates the structural evolution of the Li₅PS//Li interface at 525 K. The dashed black lines represent the initial interface positions at t = 0. Throughout the simulation, even at an elevated temperature of 1075 K, no migration of phosphorus or sulphur atoms is observed, indicating that the interfaces maintain their original structural integrity over time (unlike the γ-Li₃PS₄ and Li₇P₃S₁₁ cases). Further evidence of the high stability is provided by the positional histogram shown in Fig. 7b which reveals that after 100 ps, P and S atoms remain at their initial positions without diffusing into the lithium metal, mirroring the behaviour of Li₂S//Li system. A similar conclusion is drawn from the PDFs, as the separation distances of each pair remain almost unchanged, indicating stable atomic configurations over time (Fig. 7c).

Fig. 7 (a) Evolution of the microstructure for the Li₅PS//Li interface at 525 K. Black dashed lines indicate the initial positions of interface. (b) The positional histogram perpendicular to the interface for P and S atoms. (c) PDFs of selected atomic pairs for Li₅PS//Li at different times. (d) MSD for S atoms in Li₂S//Li, Li₇P₃S₁₁//Li, Li₅PS//Li and γ-Li₃PS₄//Li interface systems. Calculated tracer diffusivities (D*) are also displayed. The coordination numbers for (e) P atoms and (f) S atoms in the Li₅PS//Li interface structure. (g) The microstructure of the Li₅PS//Li system after 100 ps at 525 K. The local microstructures around P and S atoms are also displayed alongside bulk Li₂S and bulk Li₃P for comparison.

To quantitatively assess the interface stability, the MSD of sulphur atoms across different systems is plotted in Fig. 7d. As shown in the previous section, the instability of the γ-Li₃PS₄//Li and Li₇P₃S₁₁//Li interfaces results in the diffusion of sulphur and phosphorus atoms into the lithium metal, leading to higher MSD values and tracer diffusion coefficients. In both systems, the sulphur migration tends to decelerate after most of the sulphur atoms are consumed for SEI layer formation in a similar fashion to SEI layer growth shown in Fig. 6e. In contrast, the sulphur atoms in the Li₅PS//Li system exhibit mobility and a tracer diffusion coefficient comparable to that of the stable Li₂S//Li system indicating high stability.

The coordination numbers for phosphorus and sulphur atoms within the Li₅PS//Li system exhibit minimal variation over time (Fig. 7e and f), further supporting the conclusion that the interface remains stable. In Fig. 7g, the microstructural environment surrounding the phosphorus and sulphur atoms in the Li₅PS layer is analysed. A closer examination of individual phosphorus and sulphur atoms reveals that they are situated in Li₃P-like and Li₂S-like environments, respectively. This suggests that, unlike in the γ-Li₃PS₄ and Li₇P₃S₁₁ cases (see Fig. 6d), the phosphorus atoms in Li₅PS (and other Li₂S–Li₃P solid-solutions) are already in a 3⁻ oxidation state, which diminishes the driving force for decomposition reactions. This observation is corroborated by the energetics of the corresponding decomposition reaction of Li₅PS:

Li₅PS → Li₂S + Li₃P +0.35 eV.

The energy gain from this reaction is negligible, approximately two orders of magnitude lower than that for Li₇P₃S₁₁, further underscoring the thermodynamic stability of the Li₅PS//Li interface.

Next, we assessed the electrochemical stability of novel ternary electrolytes in contact with Li metal by computing both the decomposition electrochemical stability window (DESW) and the intrinsic electrochemical stability window (IESW), following the methodology outlined by Wagemaker et al.⁹⁷ For DESW calculations, grand potential phase diagrams were employed to evaluate the thermodynamic stability of the ternary compounds against Li metal and to determine their decomposition products as a function of applied voltage. The results, summarized in Table S4,† indicate that the DESW spans from 0 to 0.87 V, suggesting that these ternaries remain thermodynamically stable against a Li-metal anode only at low voltages. Although the stability window is relatively narrow and the decomposition into Li₂S and LiP occurs beyond 0.87 V due to delithiation, these ternaries are still promising candidates for use as protective coatings to prevent SEI formation in highly conductive electrolytes that are otherwise unstable in contact with Li metal.

The IESW, in contrast, considers the indirect decomposition of the ternaries through delithiation. To compute the IESW, we fully delithiated the ternary solid electrolytes and applied the following equation:^39,98

where V̄ represents the average voltage associated with stepwise Li removal, E_T and E_dl are the DFT-calculated total energies of the pristine and fully delithiated solid electrolytes, respectively, n_Li denotes the number of Li atoms removed, µ_Li is the chemical potential of Li, and F is the Faraday constant. The computed IESWs, as shown in Table S4,† are broader than the decomposition windows for all ternary electrolytes. This result arises from the assumption that direct decomposition is kinetically hindered, leading to indirect decomposition via delithiation. Consequently, kinetic stabilization extends the electrochemical stability window. This assumption is further corroborated by the stronger agreement between the intrinsic stability window and experimentally observed electrochemical stability in other sulphide-based solid electrolytes, such as Li₃PS₄, LGPS, Li₆PS₅Cl, Li₆PS₅Br.⁹⁷

We note that high stability in an AIMD simulation for 100 ps does not guarantee the interface stability at a much later time or in real life batteries. However, these findings indicate that the newly discovered ternaries present remarkable improvement over the well-established counterparts. For real life applications, experimental verification is required.

3.5. Ionic diffusion of Li–P–S ternaries

In the previous section, we demonstrated that novel ternary compounds exhibit stability against the lithium metal anode, due to their Li₂S-like structural framework. However, Li₂S is known for its notoriously low ionic diffusivity.^70,71 In this section, we investigate the diffusivity of Li–P–S ternary systems, with a particular emphasis on the impact of the interface on the diffusion properties.

Fig. 8 compares the AIMD lithium ion diffusion trajectories for bulk Li₂S, Li₇P₃S₁₁, and Li₅PS systems, as well as their respective interfaces with lithium metal (Li₂S//Li, Li₇P₃S₁₁//Li, and Li₅PS//Li). Li₂S systems are depicted at an elevated temperature of 1075 K, as they exhibit negligible diffusion at 525 K. In bulk Li₂S, although very limited, T–O–T type diffusion pathways are visible at 1075 K. In the Li₂S//Li interface system, on the other hand, there is an increase in diffusion, likely due to the possibility of Li atom exchange at the interface. However, the increase in diffusion is predominantly in the direction parallel to the interface, which limits its contribution to overall ionic conduction. In bulk Li₇P₃S₁₁, well-defined (yet somewhat irregular) conduction pathways are observed in agreement with previous studies,^{64–66,74,75} indicating a stable structural framework at 525 K. Conversely, in the interface structure, these pathways are disrupted due to rapid interface degradation and SEI layer formation. Notably, the formation of SEI layer impedes the movement of lithium ions which is visible as a grid-like pattern on the left side of the image for Li diffusion in Li₇P₃S₁₁//Li system. Li₅PS exhibits both high structural stability and high ionic conductivity, leading to a well-defined and similar ionic diffusion network in both bulk Li₅PS and Li₅PS//Li interface systems. The diffusion is three-dimensional and continues seamlessly (percolated) across the interface.

Fig. 8 Li⁺-ion trajectory (green) for bulk Li₂S, Li₇P₃S₁₁, and Li₅PS systems, and their respective interfaces with lithium metal (Li₂S//Li, Li₇P₃S₁₁//Li, and Li₅PS//Li). Li₇P₃S₁₁ and Li₅PS systems are presented at 525 K but Li₂S systems are depicted at 1075 K, due to negligible diffusion at 525 K. The positions of the phosphorus and sulphur atoms at t = 0 are shown in the background. Similar Li-ion trajectory plots for other systems can be found in ESI, Fig. S9.†

3.5.1. Van Hove correlation of Li⁺-ions. Van Hove correlation function is an analytical tool for probing the microscopic dynamics of ion transport in solid and liquid electrolytes. By decomposing the function into self and distinct components, ion mobility and interactions over time can be studied to understand ionic conductivity and diffusion mechanisms.

Consider a system of N particles, specifically Li atoms, with time-dependent position coordinates r_i(t), where i = 1, …, N represents the particle index, and t denotes time. The Van Hove correlation function is defined as the probability of finding a particle at position r at time t, given that a particle was located at the origin at t = 0:

where 〈 〉 denotes an ensemble average and δ is the three-dimensional Dirac delta function. G(r, t) can be decomposed into two components, conventionally termed the “self” and the “distinct” part, by distinguishing between the cases i = j and i ≠ j, respectively:

For a given distance r and time t, the self-part G_s(r, t) quantifies the likelihood that a particle has moved from its initial position by a distance r after a time interval t. Conversely, the distinct part G_d(r, t) describes the radial distribution of the remaining N − 1 particles at time t relative to the initial reference particle. Notably, G_d(r, t) reduces to the static PDF when t = 0.

In this study, we calculated the distinct part of the Van Hove correlation function for the charge carrier Li⁺ ions. If the system possesses well-defined lattice sites for Li atoms and some of these sites remain occupied throughout the AIMD simulation, pronounced peaks emerge in the function. However, if the initial lattice structure is disrupted and the well-defined sites are lost, the reference positions of Li atoms become less meaningful, leading to a broadening and eventual disappearance of the peaks.

Fig. 9 presents the distinct part of the Van Hove correlation function for Li₂S, Li₇P₃S₁₁, and Li₅PS systems. For bulk Li₂S at 525 K, the lack of lithium diffusion results in a plot that closely resembles the PDF over time, exhibiting the first and second coordination shells at approximately 3 Å and 4 Å, respectively. However, at the elevated temperature of 1075 K, the lithium atoms in Li₂S acquire sufficient energy to diffuse, leading to a pronounced peak in G_d(r, t) near r = 0, indicating highly correlated Li⁺-ion motions and a high probability that a vacated Li site is quickly occupied by another Li atom. A similar behaviour is observed in bulk Li₇P₃S₁₁, bulk Li₅PS, and Li₅PS//Li systems at 525 K. In Li₇P₃S₁₁//Li system, however, the peak near r = 0 is weak and diffuse. This can be attributed to degradation at the interface and the formation of the SEI layer, which disrupts the well-defined initial Li occupation sites. Additionally, the coordination shells in both bulk Li₇P₃S₁₁ and Li₇P₃S₁₁//Li systems appear less distinct. In contrast, both bulk Li₅PS and Li₅PS//Li systems exhibit pronounced peaks near r = 0 and well-defined coordination shells due to high stability and diffusivity as discussed in previous sections.

Fig. 9 The distinct part of the Van Hove correlation function for Li atoms in (a) bulk Li₂S at 525 K, (b) bulk Li₂S at 1075 K, (c) bulk Li₇P₃S₁₁ at 525 K, (d) Li₇P₃S₁₁//Li at 525 K, (e) Li₅PS at 525 K and, (f) Li₅PS//Li at 525 K. Similar Van Hove correlation plots for other bulk and interface systems can be found in ESI, Fig. S10.†

3.5.2. Diffusion at the interface. To investigate the nanoscale diffusion properties of electrolytes, particularly at their interfaces with Li metal, we developed a Python code that tracks the number of Li⁺-ions crossing an imaginary zone, as illustrated in Fig. 10a, in both directions along the c-axis. The c-axis in the bulk systems are aligned with specific Miller indices to ensure the same orientation as in the interface structure. The zone thickness is set to 3 Å to avoid counting oscillatory motions of Li⁺-ions at a single site, which do not represent actual site-to-site jumps. To examine the impact of interfacing the electrolyte with lithium metal anode, the zone is positioned within the bulk electrolyte (referred to as “bulk” zone), at the interface (referred to as “interface” zone), and at the centre (referred to as “centre” zone) of the electrolyte in the interface system. We note that although these diffusion events are sampled under fixed-temperature conditions, as some of the interfaces degrade during the simulations, not all of them necessarily reach thermodynamic equilibrium or a steady state. Furthermore, no external electric field is applied to the ions in these AIMD simulations, unlike in real-life charging/discharging processes.

Fig. 10 Local diffusion calculations for Li⁺-ions. (a) Bulk (black), interface (red) and centre (blue) imaginary zones that are used to calculate the flux of ions are depicted in bulk Li₅PS and Li₅PS//Li interface structures. (b) Bidirectional ion flux at 525 K through bulk, interface and centre zones for bulk Li₅PS and Li₇P₃S₁₁ along with Li₅PS//Li and Li₇P₃S₁₁//Li interface systems. (c) Directional flux at 525 K into and out of the electrolyte in Li₅PS//Li and Li₇P₃S₁₁//Li interface models.

In Fig. 10b, we compare the bidirectional flux of Li⁺-ions—the total count of Li atoms passing through the zone in either direction—for bulk Li₅PS and bulk Li₇P₃S₁₁, as well as their interfaces with the lithium metal anode at 525 K. For Li₅PS, the flux at both centre and interface zones is very similar to that observed at bulk zone, with a constant increase over time. This indicates the absence of significant internal resistance to Li⁺-ion flux at the interface. In contrast, for Li₇P₃S₁₁, the flux at the “interface” is lower than that at the “centre”, suggesting the presence of resistance, likely due to the formation of the SEI layer that impedes Li⁺-ion flux. However, the fluxes at both interface and centre zones are significantly higher than at bulk zone, indicating the presence of a different diffusion mechanism in the Li₇P₃S₁₁//Li system than its bulk counterpart.

Fig. 10c presents the directional flux of Li⁺-ions—the total count of Li atoms passing through the zone in one direction—into and out of the electrolytes in Li₅PS//Li and Li₇P₃S₁₁//Li systems at 525 K. Initially, the slope of the Li⁺-ion flux into the electrolyte is much steeper than the flux out, indicating a rapid influx of Li atoms into the electrolyte in Li₇P₃S₁₁//Li system. After 40 ps, the slopes of the fluxes in each direction become similar, suggesting a steady state. The initial surge of Li ions into the electrolyte is attributed to the porous nature of the material. In our simulations, in bulk Li₇P₃S₁₁, ionic conduction occurs via Li ion jumps from one site to another. However, when Li₇P₃S₁₁ is interfaced with the Li metal, the excess Li atoms rapidly infiltrate the pores of Li₇P₃S₁₁, not only leading to the formation of an SEI layer but also providing an alternative conduction pathway through vacant interstitial sites, leading to a change in the principle conduction mechanism. In contrast, the structural framework of Li₅PS is robust and without such pores, which dictates that the conduction mainly occurs through stable vacant interstitial sites, as discussed in Section 3.1. Consequently, the flux of Li ions in each direction through the interface is balanced and the steady-state is observed throughout the simulation.

4. Conclusions

Our comprehensive DFT simulations have revealed the diffusion mechanisms and the underlying reasons for the high ionic conductivity in novel Li–P–S ternary electrolytes, Li₇PS₂, Li₅PS, Li₈P₂S, and Li₁₁P₃S. These materials, exhibiting ionic conductivities on par with the well-known highly-conductive Li₇P₃S₁₁, represent a promising class of electrolytes for solid-state batteries. Among the newly studied ternaries, Li₅PS stands out with the highest conductivity, a property we attribute to the homogenous distribution of phosphorus and sulphur atoms within its structure (due to the 1 : 1 P to S ratio). This balanced atomic arrangement facilitates more efficient lithium-ion transport, reducing energy barriers and enabling faster ionic movement through the lattice.

Our investigation also extended to the thermodynamic stability of these ternaries in contact with lithium metal anodes, a crucial factor for practical Li-metal battery applications. Compared to Li₇P₃S₁₁, which rapidly forms an SEI layer due to the reduction of phosphorus from an oxidation state of 5⁺ to 3⁻, the novel Li–P–S ternaries demonstrate remarkable stability. The SEI formation rate in Li₇P₃S₁₁ is a time-dependent process, where the thickness of the interphase layer grows proportionally to the square root of time, progressively increasing the internal resistance and reducing the battery's overall efficiency. The stability of novel ternaries is rooted in their unique Li₂S-like structural framework, which inherently positions phosphorus atoms in the 3⁻ oxidation state, thereby eliminating the main (chemical) driving force for SEI formation. The absence of SEI layer is particularly significant, as it means that no additional resistance is introduced at the electrolyte–anode interface, which is a common issue in many solid-state battery systems. These finding not only highlights the suitability of novel ternaries as high-performance SEs but also underscore the importance of structural design in developing next-generation battery materials. By avoiding the issues associated with SEI formation, the novel Li–P–S ternaries can maintain high ionic conductivity and stability over extended periods, making them ideal candidates as electrolytes and coating materials for use in solid-state lithium batteries.

In summary, our study provides a detailed understanding of the physical mechanisms driving the superior performance of these new Li–P–S ternary electrolytes. The insights gained from our work pave the way for the development of advanced solid-state batteries that combine high ionic conductivity with robust chemical and electrochemical stability. These findings suggest that the novel ternaries, particularly Li₅PS, could play a crucial role in the future of energy storage technologies, offering a pathway to safer, more efficient, and longer-lasting batteries.

Data availability

DFT-optimised structure files of bulk, interface and surfaces (with minimum surface formation energies) are available in the GitHub repository: https://www.github.com/bkarasulu/LiPS-interfaces-data.git.

Author contributions

H. S. S. contributed in conceptualization, data curation, formal analysis, investigation, validation, visualization, methodology, writing, reviewing and editing of the draft. B. K. contributed in conceptualization, software, resources, funding acquisition, supervision, project administration, methodology, writing, reviewing and editing of the draft.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

B. K. and H. S. S. acknowledge support from the UK Research and Innovation (UKRI)'s Engineering and Physical Sciences Research Council (EPSRC) through the Early Career Fellowship (EP/T026138/1). Computing facilities were provided by the Scientific Computing Research Technology Platform (SC-RTP) at the University of Warwick. DFT calculations were performed using the Avon and Sulis HPC platforms, and Sulis is funded by the EPSRC Grant EP/T022108/1 and the HPC Midlands + consortium. This work also used the ARCHER2 UK National Supercomputing Service (https://www.archer2.ac.uk) for performing the DFT calculations (through the Project e848).

Notes and references

1. X. Feng, D. Ren, X. He and M. Ouyang, Joule, 2020, 4, 743–770.
2. Z. Chen, R. Xiong, J. Lu and X. Li, Appl. Energy, 2018, 213, 375–383.
3. D. P. Finegan, M. Scheel, J. B. Robinson, B. Tjaden, I. Hunt, T. J. Mason, J. Millichamp, M. Di Michiel, G. J. Offer and G. Hinds, Nat. Commun., 2015, 6, 6924.
4. T. Krauskopf, F. H. Richter, W. G. Zeier and J. Janek, Chem. Rev., 2020, 120, 7745–7794.
5. X.-B. Cheng, R. Zhang, C.-Z. Zhao and Q. Zhang, Chem. Rev., 2017, 117, 10403–10473.
6. C. Fang, B. Lu, G. Pawar, M. Zhang, D. Cheng, S. Chen, M. Ceja, J. M. Doux, H. Musrock, M. Cai, B. Liaw and Y. S. Meng, Nat. Energy, 2021, 6, 987–994.
7. X. Liang, Q. Pang, I. R. Kochetkov, M. S. Sempere, H. Huang, X. Sun and L. F. Nazar, Nat. Energy, 2017, 2, 1–7.
8. D. Lin, Y. Liu and Y. Cui, Nat. Nanotechnol., 2017, 12, 194–206.
9. P. Xue, S. Liu, X. Shi, C. Sun, C. Lai, Y. Zhou, D. Sui, Y. Chen and J. Liang, Adv. Mater., 2018, 30, 1804165.
10. X.-B. Cheng, R. Zhang, C.-Z. Zhao, F. Wei, J.-G. Zhang and Q. Zhang, Adv. Sci., 2016, 3, 1500213.
11. W. Xu, J. Wang, F. Ding, X. Chen, E. Nasybulin, Y. Zhang and J. G. Zhang, Energy Environ. Sci., 2014, 7, 513–537.
12. P. Albertus, S. Babinec, S. Litzelman and A. Newman, Nat. Energy, 2018, 3, 16–21.
13. N. Ohta, K. Takada, L. Zhang, R. Ma, M. Osada and T. Sasaki, Adv. Mater., 2006, 18, 2226–2229.
14. K. Takada, T. Inada, A. Kajiyama, H. Sasaki, S. Kondo, M. Watanabe, M. Murayama and R. Kanno, Solid State Ionics, 2003, 158, 269–274.
15. L. Lin, F. Liang, K. Zhang, H. Mao, J. Yang and Y. Qian, J. Mater. Chem. A, 2018, 6, 15859–15867.
16. F. Lalère, J.-B. Leriche, M. Courty, S. Boulineau, V. Viallet, C. Masquelier and V. Seznec, J. Power Sources, 2014, 247, 975–980.
17. H. Kitaura and H. Zhou, Sci. Rep., 2015, 5, 13271.
18. J. Janek and W. G. Zeier, Nat. Energy, 2016, 1, 1–4.
19. R. Murugan, V. Thangadurai and W. Weppner, Angew. Chem., Int. Ed., 2007, 46, 7778–7781.
20. P. Knauth, Solid State Ionics, 2009, 180, 911–916.
21. J. C. Bachman, S. Muy, A. Grimaud, H.-H. Chang, N. Pour, S. F. Lux, O. Paschos, F. Maglia, S. Lupart and P. Lamp, Chem. Rev., 2016, 116, 140–162.
22. C. Dietrich, D. A. Weber, S. J. Sedlmaier, S. Indris, S. P. Culver, D. Walter, J. Janek and W. G. Zeier, J. Mater. Chem. A, 2017, 5, 18111–18119.
23. S. S. Berbano, M. Mirsaneh, M. T. Lanagan and C. A. Randall, Int. J. Appl. Glass Sci., 2013, 4, 414–425.
24. 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.
25. S. V. Patel, S. Banerjee, H. Liu, P. Wang, P. H. Chien, X. Feng, J. Liu, S. P. Ong and Y. Y. Hu, Chem. Mater., 2021, 33, 1435–1443.
26. Y. Kato, S. Hori, T. Saito, K. Suzuki, M. Hirayama, A. Mitsui, M. Yonemura, H. Iba and R. Kanno, Nat. Energy, 2016, 1, 16030.
27. S. Strangmüller, H. Eickhoff, D. Müller, W. Klein, G. Raudaschl-Sieber, H. Kirchhain, C. Sedlmeier, V. Baran, A. Senyshyn, V. L. Deringer, L. Van Wüllen, H. A. Gasteiger and T. F. Fässler, J. Am. Chem. Soc., 2019, 141, 14200–14209.
28. T. M. F. Restle, C. Sedlmeier, H. Kirchhain, W. Klein, G. Raudaschl‐Sieber, V. L. Deringer, L. van Wüllen, H. A. Gasteiger and T. F. Fässler, Angew. Chem., 2020, 132, 5714–5723.
29. S. Strangmüller, H. Eickhoff, G. Raudaschl-Sieber, H. Kirchhain, C. Sedlmeier, L. Van Wüllen, H. A. Gasteiger and T. F. Fässler, Chem. Mater., 2020, 32, 6925–6934.
30. T. M. F. Restle, C. Sedlmeier, H. Kirchhain, W. Klein, G. Raudaschl-Sieber, L. Van Wüllen and T. F. Fässler, Chem. Mater., 2021, 33, 2957–2966.
31. L. Toffoletti, H. Kirchhain, J. Landesfeind, W. Klein, L. van Wüllen, H. A. Gasteiger and T. F. Fässler, Chem.—Eur. J., 2016, 22, 17635–17645.
32. H. Eickhoff, S. Strangmüller, W. Klein, H. Kirchhain, C. Dietrich, W. G. Zeier, L. van Wüllen and T. F. Fässler, Chem. Mater., 2018, 30, 6440–6448.
33. S. Strangmüller, H. Eickhoff, W. Klein, G. Raudaschl-Sieber, H. Kirchhain, T. Kutsch, V. Baran, A. Senyshyn, L. van Wüllen, H. A. Gasteiger and T. F. Fässler, J. Mater. Chem. A, 2021, 9, 15254–15268.
34. W. D. Richards, L. J. Miara, Y. Wang, J. C. Kim and G. Ceder, Chem. Mater., 2016, 28, 266–273.
35. F. Han, Y. Zhu, X. He, Y. Mo and C. Wang, Adv. Energy Mater., 2016, 6, 1501590.
36. S. P. Ong, Y. Mo, W. D. Richards, L. Miara, H. S. Lee and G. Ceder, Energy Environ. Sci., 2013, 6, 148–156.
37. Y. Zhu, X. He and Y. Mo, Advanced Science, 2017, 4, 1600517.
38. S. Wenzel, S. Randau, T. Leichtweiß, D. A. Weber, J. Sann, W. G. Zeier and J. Janek, Chem. Mater., 2016, 28, 2400–2407.
39. A. F. Harper, M. L. Evans, J. P. Darby, B. Karasulu, C. P. Koçer, J. R. Nelson and A. J. Morris, Johnson Matthey Technol. Rev., 2020, 64, 103–118.
40. Y. Wang, W. D. Richards, S. P. Ong, L. J. Miara, J. C. Kim, Y. Mo and G. Ceder, Nat. Mater., 2015, 14, 1026–1031.
41. A. Van Der Ven, Z. Deng, S. Banerjee and S. P. Ong, Chem. Rev., 2020, 120, 6977–7019.
42. M. Pasta, D. Armstrong, Z. L. Brown, J. Bu, M. R. Castell, P. Chen, A. Cocks, S. A. Corr, E. J. Cussen, E. Darnbrough, V. Deshpande, C. Doerrer, M. S. Dyer, H. El-Shinawi, N. Fleck, P. Grant, G. L. Gregory, C. Grovenor, L. J. Hardwick, J. T. S. Irvine, H. J. Lee, G. Li, E. Liberti, I. McClelland, C. Monroe, P. D. Nellist, P. R. Shearing, E. Shoko, W. Song, D. S. Jolly, C. I. Thomas, S. J. Turrell, M. Vestli, C. K. Williams, Y. Zhou and P. G. Bruce, JPhys Energy, 2020, 2, 032008.
43. A. R. Oganov, C. J. Pickard, Q. Zhu and R. J. Needs, Nat. Rev. Mater., 2019, 4, 331–348.
44. Z. Zhu, I. H. Chu and S. P. Ong, Chem. Mater., 2017, 29, 2474–2484.
45. R. P. Rao, H. Chen and S. Adams, Chem. Mater., 2019, 31, 8649–8662.
46. 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.
47. N. L. Marana, S. Casassa, M. F. Sgroi, L. Maschio, F. Silveri, M. D'Amore and A. M. Ferrari, Langmuir, 2023, 39, 18797–18806.
48. N. D. Lepley and N. A. W. Holzwarth, Phys. Rev. B:Condens. Matter Mater. Phys., 2015, 92, 214201.
49. H. Chen, A. Pei, D. Lin, J. Xie, A. Yang, J. Xu, K. Lin, J. Wang, H. Wang, F. Shi, D. Boyle and Y. Cui, Adv. Energy Mater., 2019, 9, 1900858.
50. INTERFACER: A Python-script for automated generation of atomistic models for solid-solid interfaces., Freely available under GNU license at https://github.com/bkarasulu/INTERFACER-code .
51. G. Kresse and J. Hafner, Phys. Rev. B:Condens. Matter Mater. Phys., 1993, 47, 558.
52. G. Kresse and J. Furthmüller, Comput. Mater. Sci., 1996, 6, 15–50.
53. G. Kresse and J. Furthmüller, Phys. Rev. B:Condens. Matter Mater. Phys., 1996, 54, 11169.
54. P. E. Blöchl, Phys. Rev. B:Condens. Matter Mater. Phys., 1994, 50, 17953.
55. G. Kresse and D. Joubert, Phys. Rev. B:Condens. Matter Mater. Phys., 1999, 59, 1758.
56. J. P. Perdew and W. Yue, Phys. Rev. B:Condens. Matter Mater. Phys., 1986, 33, 8800.
57. K. Momma and F. Izumi, J. Appl. Crystallogr., 2011, 44, 1272–1276.
58. A. Stukowski, Modell. Simul. Mater. Sci. Eng., 2009, 18, 15012.
59. J. Ding, M. K. Gupta, C. Rosenbach, H. M. Lin, N. C. Osti, D. L. Abernathy, W. G. Zeier and O. Delaire, Nat. Phys., 2025, 21, 118–125.
60. T. Krauskopf, C. Pompe, M. A. Kraft and W. G. Zeier, Chem. Mater., 2017, 29, 8859–8869.
61. T. Krauskopf, S. Muy, S. P. Culver, S. Ohno, O. Delaire, Y. Shao-Horn and W. G. Zeier, J. Am. Chem. Soc., 2018, 140, 14464–14473.
62. A. Kuhn, J. Köhler and B. V. Lotsch, Phys. Chem. Chem. Phys., 2013, 15, 11620–11622.
63. K. Homma, M. Yonemura, T. Kobayashi, M. Nagao, M. Hirayama and R. Kanno, Solid State Ionics, 2011, 182, 53–58.
64. Y. Onodera, K. Mori, T. Otomo, M. Sugiyama and T. Fukunaga, J. Phys. Soc. Jpn., 2012, 81, 044802.
65. D. Chang, K. Oh, S. J. Kim and K. Kang, Chem. Mater., 2018, 30, 8764–8770.
66. H. Yamane, M. Shibata, Y. Shimane, T. Junke, Y. Seino, S. Adams, K. Minami, A. Hayashi and M. Tatsumisago, Solid State Ionics, 2007, 178, 1163–1167.
67. R. J. Friauf, J. Appl. Phys., 1962, 33, 494–505.
68. X. He, Y. Zhu, A. Epstein and Y. Mo, npj Comput. Mater., 2018, 4, 18.
69. M. K. Gupta, R. Mittal, B. Singh, O. Delaire, S. N. Achary, S. Rols, A. K. Tyagi and S. L. Chaplot, Phys. Rev. B, 2021, 103, 174109.
70. S. P. Jand, Q. Zhang and P. Kaghazchi, Sci. Rep., 2017, 7, 5873.
71. F. Altorfer, W. Buhrer, I. Anderson, O. Scharpf, H. Bill, P. L. Carron and H. G. Smith, Phys. B: Condens. Matter, 1992, 180–181, 795–797.
72. W. Buehrer, F. Altorfer, J. Mesot, H. Bill, P. Carron and H. G. Smith, J. Phys.:Condens. Matter, 1991, 1055–1064.
73. X. He, Y. Zhu and Y. Mo, Nat. Commun., 2017, 8, 15893.
74. Y. Onodera, K. Mori, T. Otomo, A. C. Hannon, S. Kohara, K. Itoh, M. Sugiyama and T. Fukunaga, J. Phys. Soc. Jpn., 2010, 79, 87–89.
75. K. Xiong, R. C. Longo, S. Kc, W. Wang and K. Cho, Comput. Mater. Sci., 2014, 90, 44–49.
76. J. S. Park, Y. K. Jung, K. T. Butler and A. Walsh, J. Phys.: Energy, 2019, 1, 016001.
77. K. T. Butler, G. Sai Gautam and P. Canepa, npj Comput. Mater., 2019, 5, 19.
78. W. Xie, Z. Deng, Z. Liu, T. Famprikis, K. T. Butler and P. Canepa, Adv. Energy Mater., 2024, 14, 2304230.
79. P. Canepa, J. A. Dawson, G. Sai Gautam, J. M. Statham, S. C. Parker and M. S. Islam, Chem. Mater., 2018, 30, 3019–3027.
80. B. Gao, R. Jalem and Y. Tateyama, ACS Appl. Mater. Interfaces, 2020, 12, 16350–16358.
81. T. Jia, D. J. Senor and Y. Duan, J. Nucl. Mater., 2021, 555, 153111.
82. H. Wu, J. Zhang, Y. Wang, J. Shang and Y. Jiang, Coatings, 2023, 13, 1203.
83. B. Karasulu, S. P. Emge, M. F. Groh, C. P. Grey and A. J. Morris, J. Am. Chem. Soc., 2020, 142, 3132–3148.
84. J. Haruyama, K. Sodeyama, L. Han, K. Takada and Y. Tateyama, Chem. Mater., 2014, 26, 4248–4255.
85. J. Wang, A. A. Panchal, G. Sai Gautam and P. Canepa, J. Mater. Chem. A, 2022, 10, 19732–19742.
86. R. Tran, Z. Xu, B. Radhakrishnan, D. Winston, W. Sun, K. A. Persson and S. P. Ong, Sci. Data, 2016, 3, 160080.
87. I. D. Seymour and A. Aguadero, J. Mater. Chem. A, 2021, 9, 19901–19913.
88. N. D. Lepley, N. A. W. Holzwarth and Y. A. Du, Phys. Rev. B:Condens. Matter Mater. Phys., 2013, 88, 104103.
89. I. H. Chu, H. Nguyen, S. Hy, Y. C. Lin, Z. Wang, Z. Xu, Z. Deng, Y. S. Meng and S. P. Ong, ACS Appl. Mater. Interfaces, 2016, 8, 7843–7853.
90. T. L. Einstein, in Handbook of Crystal Growth, ed. T. Nishinaga, Elsevier, 2 edn, 2015, pp. 215–264.
91. S. Miracle-Sole, in Mathematical Results in Statistical Mechanics, Satellite Colloquium of STATPHYS, 1999, vol. 20, pp. 83–101.
92. W. D. Nix and H. Gao, Scr. Mater., 1998, 39, 1653–1661.
93. C.-E. Kim, Y.-J. Tak, K. T. Butler, A. Walsh and A. Soon, Phys. Rev. B:Condens. Matter Mater. Phys., 2015, 91, 85307.
94. L. E. Camacho-Forero and P. B. Balbuena, J. Power Sources, 2018, 396, 782–790.
95. T. Yamada, S. Ito, R. Omoda, T. Watanabe, Y. Aihara, M. Agostini, U. Ulissi, J. Hassoun and B. Scrosati, J. Electrochem. Soc., 2015, 162, A646–A651.
96. N. A. W. Holzwarth, N. D. Lepley and Y. A. Du, J. Power Sources, 2011, 196, 6870–6876.
97. T. K. Schwietert, A. Vasileiadis and M. Wagemaker, JACS Au, 2021, 1, 1488–1496.
98. A. Urban, D.-H. Seo and G. Ceder, npj Comput. Mater., 2016, 2, 16002.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d5ta00585j

---