Atomistic insights into the chemical stability and ionic transport at Li-metal/Li-argyrodite interfaces

Abstract

All solid-state batteries (ASSBs) based on solid-state electrolytes (SSEs) are a novel Li-ion battery technology with the potential of enhanced safety, longer lifetimes, and increased energy density when coupled with the Li-metal anode. Li-argyrodite (Li₆PS₅Cl) is a promising SSE with high ionic conductivity, produced using cheap and sustainable precursors, and therefore of interest to both academia and industry. Like many other sulfide-based SSEs, it is however unstable against Li-metal. Using ab initio and machine-learning methods, we simulate three representative Li-metal/Li-argyrodite interface models to investigate whether the exact surface termination affects the chemical stability and ion transport capability. We present a systematic approach to create low-energy interfaces by screening 28 low Miller-index surface terminations of Li-argyrodite and coupling them with Li-metal. Custom-made machine-learned interatomic potentials trained on ab initio data enable the simulation of large interface models with over 2000 atoms for 5 ns. We find that all three interfaces decompose into an amorphous solid-electrolyte interphase (SEI) layer, consisting of Li₃P, Li₂S and LiCl, which then crystallises into an antifluorite phase Li₂S_1−x−yP_xCl_y; {x, y = 0.14–0.15}. Calculating the flux of Li-ions across a fixed membrane, we show that the crystalline SEI layer is a poor ion conductor, similar to Li₂S. While all three interfaces form the same crystalline SEI layer, the exact rates of the decomposition and crystallisation depend on the actual surface composition. These atomic-level insights will be practical to control the SEI formation in sulphide-based SSEs and others.

---

1 Introduction

All-solid-state batteries (ASSBs) are a new battery technology that has several benefits over conventional liquid Li-ion batteries (LLIBs). Replacing the flammable liquid electrolyte with a solid-state electrolyte (SSE) alleviates some of the safety concerns surrounding LLIBs and can lead to longer lifetimes as well as higher energy and power densities. Moreover, ASSBs might enable the use of the Li-metal anode, which would further increase the energy density due to its high theoretical capacity of 3860 mA h g⁻¹ and low electrochemical potential of −3.04 V versus the standard hydrogen electrode.^1–4

There are, however, several challenges that still need to be overcome before commercialisation, most of which are related to the interface between solid-state components. One problem is the chemical and electrochemical stability of the interfaces, as many SSE materials are easily reduced at low potentials, which leads to solid electrolyte interface (SEI) formation and degradation of the electrolyte. Another important issue is mechanical failure, such as the formation and propagation of cracks, and delamination of the interface, which is mainly caused by the volume expansion of the electrode and hinders Li-ion transport, thus affecting the performance and lifetime of ASSBs.^1–3,5,6 To find suitable solutions to these problems, a good understanding of the underlying structure of solid–solid interfaces, the chemical stability, and the pathway of diffusion across the interface is required.¹

There is currently a strong interest in sulfide SSEs, with several companies attempting to commercialise sulfide-containing ASSBs.⁷ Amongst the sulfide SSEs, the Li-argyrodite family, which was first reported by Deiseroth et al. in 2008 and has a composition of Li₆PS₅X (X = Cl, Br, I),⁸ has many desirable properties. The chloride equivalent Li₆PS₅Cl has a high ionic conductivity of 10⁻³ to 10⁻² S cm⁻¹ and can be synthesised using cheaper precursors than other SSEs of comparable conductivity, such as Li₁₀GeP₂S₁₂.⁹ Due to the aforementioned benefits of the Li-metal anode, the Li-metal/Li-argyrodite interface has gained significant attention in recent years both in academia and industry.

Computational methods are well suited to study bulk material properties as well as solid–solid interface phenomena and gain insight at the atomistic level. The bulk structure of Li-argyrodite has been extensively investigated computationally.^10–16 A few early computational studies investigated the Li-metal/Li-argyrodite interfaces, using first-principles methods such as Density Functional Theory (DFT) and Ab initio Molecular Dynamics (AIMD). Golov et al.¹⁷ determined the migration barriers of four different Li-metal/Li-argyrodite interfaces and Cheng et al.¹⁸ investigated the decomposition reactions occurring at the interface. However, the high computational cost of ab initio methods limits the size of the simulation models to a few 100 atoms and simulation times to below nanoseconds, which is not enough to study dynamic processes at solid–solid interfaces. With Machine Learning Interatomic Potentials (MLIPs) becoming mature enough to handle complex, multi-species systems, the time- and length scales required are now accessible while keeping the accuracy of DFT. Wang et al.¹⁹ used MLIPs to study the ion transport at interfaces between Li-metal and the reaction products of Li-argyrodite. During the preparation of this manuscript, two new studies have been published by Chaney et al.²⁰ and Ding et al.²¹ using MLIPs to drive long MD simulations (up to 35 ns for a 21 120 atom cell) and study larger systems (up to three-million-atom cells for 1 ns) of the Li-metal/Li-argyrodite interface. This allowed them to simulate the early stages of SEI formation and observe the transition from an amorphous to a new crystalline disordered phase (Li₂S_0.72P_0.14Cl_0.14). Both papers focused on the stability and SEI formation, and while Ding et al.²¹ calculated the diffusion coefficient of bulk Li₂S_0.72P_0.14Cl_0.14, they did not investigate the ion transport across the interface and how it is affected by the formation of the SEI. Furthermore, both groups investigated only one interface system (Li(110)/Li₅PS₅Cl(110) by Chaney et al. and Li(110)/Li₅PS₅Cl(100) by Ding et al.), without giving a specific reason for their choice.

This paper introduces a systematic computational approach to create representative atomistic models of solid–solid interfaces. This is done by finding the energetically most favourable Li-argyrodite surfaces and using our in-house code INTERFACER to generate interface models considering various Miller planes of Li-metal that give the lowest mismatch between the two lattices. Three interface structures were created from the three surfaces with the lowest formation energies to investigate whether the choice of surface termination has an effect on the SEI formation and the diffusion mechanism across the interfaces. Three MACE MLIP models for each of the interfaces were trained using the bulk structures of both electrode and electrolyte, and the respective interface structures, which were then used to run MD simulations of cells with over 2000 atoms for 5 ns, with no external driving forces applied. We discover that for all three interfaces, the electrolyte decomposes into an amorphous layer of Li₃P, Li₂S, and LiCl, before crystallizing into an antifluorite structure with a composition of Li₂S_1−x−yP_xCl_y; {x, y = 0.14–0.15}. However, the exact rates of the decomposition and crystallization vary for the different interfaces investigated. We also find that the Li-ion flux decreases by two orders of magnitude during the SEI formation and drops to 0 in the crystalline phase, showing that this is the bottleneck in the diffusion after the crystalline SEI has formed, and main reason for the high interfacial resistance generally observed for sulphide-based SSEs.

2 Computational methods

A model of the disordered bulk Li-argyrodite structure was created from the ordered structure. Configuration enumeration analysis was performed on the unit cell of the ordered structure to find the energetically most favourable arrangement of S and Cl atoms in the disordered system. Exchanging two of the 4 S atoms on the 4c site with Cl and then exchanging 2 Cl atoms with S results in two symmetrically inequivalent structures. Performing a geometry optimisation on both structures we chose the one with the lower energy to take forward. To validate it accurately captures the diffusion properties AIMD simulations were run on the ordered and disordered structure at 600 K. 28 surface cuts of Li-argyrodite were created from the disordered bulk model and relaxed using DFT. The three energetically most favourable surfaces were then interfaced with Li-metal, using our in-house code INTERFACER. The code systematically displaces and rotates the two surfaces to find the alignment that gives the highest interfacial area and the lowest lattice mismatch. Both surfaces are then put together into a supercell with the a and b parameters being an average of the two individual surface cell parameters. The resulting interface structures were relaxed with DFT. Some preliminary AIMD simulations were run for all three interfaces at 400 K and 1000 K, which were later also used to create a test dataset of 500 structures for each interface. For the MLIP training dataset, AIMD simulations were run for the 3 × 3 × 3 supercell of bcc Li-metal, as well as of the decomposition products (Li₃P, Li₂S, and LiCl). All other training structures were obtained from our active learning procedure, detailed in Section 2.2. The ordered Li-argyrodite structure, bcc Li-metal, as well as the decomposition structures, were taken from the Materials Project (Li₆PS₅Cl: mp-985592, Li: mp-135, Li₃P: mp-736, Li₂S: mp-1153, and LiCl: 22905).²² The specific settings for the first-principle calculations and the training of the MLIP models are given below.

2.1 First-principles – DFT and AIMD

DFT and AIMD calculations were carried out using the Vienna Ab initio Simulation Package (version 6.4)^23–25 with the Perdew–Burke–Ernzerhof (PBE)²⁶ exchange–correlation function. We employed the Projector-Augmented-Wave formalism,²⁷ with only the Li(s¹), P(s²p³), S(s²p⁴), and Cl(s²p⁵) electrons being treated explicitly. The cutoff energy for the DFT calculations (geometry optimisations and single point energy calculations) was set to 520 eV, which was dropped to 400 eV for the AIMD simulations. The Γ-point centred k-point grid was generated automatically with a spacing of 0.05 × 2π = 0.32 Å⁻¹. During geometry optimistion, the cell parameters and atomic positions were fully relaxed using the conjugate gradient algorithm until the change in forces reached 0.01 eV Å⁻¹. The AIMD simulations were performed in the NVT ensemble using a Nose–Hoover thermostat. A Nose mass corresponding to a period of 40 time steps (SMASS = 0) was employed. A time-step of 2 fs was used for all simulations except for the three interface structures, where a 1 fs time-step had to be used to ensure energy conservation. For the smaller bulk Li-argyrodite cell (52 atoms) and Li-metal cell (54 atoms), we used the automatically generated k-point grid, but otherwise the AIMD simulations were run with Γ-point only. Charge partitioning of the interface structures was performed using the DDEC6 code.²⁸ The Radial Distribution Functions (RDFs) were calculated with a 4 Å cutoff, while the neighbour lists for the calculations of the coordination numbers, SEI thickness, and crystallisation were created using a cutoff of 3 Å.

2.2 MLIP training

MACE,^29,30 a message passing neural network, was used for the construction of the MLIPs. To create a representative dataset, we used our in-house active-learning algorithm, STING³¹ shown in Fig. 1.

Fig. 1 Workflow of the active learning algorithm.

An initial dataset consisting of Li-argyrodite bulk, Li-metal bulk and interface structures taken from the DFT relaxation was provided. Here, we passed 300 structures from AIMD simulation of the bulk structure at 400 K, 600 K, 800 K, and 1000 K (75 from each temperature), 100 structures from an AIMD simulation of Li-bulk with temperature ramping between 200 K and 1200 K, and 50 structures from the geometry optimisation of the respective interface (using the same DFT settings as the AIMD runs). As starting structure for the MD runs, we use the final structure from the geometry optimisation of the respective interface. An initial MD run was performed at 400 K, 600 K, 800 K, and 1000 K for 4000 steps using the MACE foundation-model OMAT,³² based on Meta's Open Materials 2024 (OMat24) dataset³³ from which 201 structures as well as a validation set of 31 structures are taken equidistant and their energies and forces evaluated with DFT (VASP). The MD was run via the Atomic Simulation Environment (ASE)³⁴ in the NVT ensemble using the Langevin thermostat with a timestep of 1 fs and a friction coefficient of 0.1 fs⁻¹. For the DFT single-point calculations the same settings were used as for the AIMD simulations (using Γ-point only) to keep all the DFT energies and forces in the dataset consistent. The updated dataset is used to fine-tune a committee of three MACE models on the OMAT dataset, using MACE freeze.³⁵ The freeze parameter was set to 6, with energy, forces, and stress weights being 10.0, 1000.0, and 100.0. We employ a learning rate of 0.005, running 300 epochs with a batch size of 16, and enabling the exponential moving average (EMA). Stage two optimisation with revised weights of 1000.0, 100.0, and 100.0 for energy, forces and stress, respectively, is started after 200 epochs. All other parameters are either set to their default value, or determined by fine-tuning (i.e. universal loss, a radial cutoff radius of 6 Å, and 128 equivariant messages). The isolated atoms energies (E₀) were calculated by running spin-polarised single-point DFT calculations with VASP for each species, with the settings otherwise being the same as the AIMD settings that are used for the dataset structures.

Three MD simulations are run for each model for 6000 steps with the same NVT settings as detailed above, to determine the committee uncertainty and, if this is not yet satisfactory, select the 200 structures with the highest uncertainty in their force components. These are added to the dataset, and the steps (1) fine-tune a committee of three MACE models, (2) run committee MD and determine uncertainty, and (3) select 200 structures with the highest uncertainty, are repeated until the uncertainty in the forces is converged (see Fig. S1), which in our case was after 6 iterations. A total of 1401 structures from the initial MD run and the 6 iterations were added to the initial dataset of 450 structures, resulting in a dataset of 1851 structures.

To ensure the MD is stable, even for larger systems, we added structures of AIMD runs of the decomposition products, Li₃P, Li₂S, and LiCl. Each AIMD simulation was run with a temperature ramping from 200 K to 1200 K, and 100 structures were selected equidistant for each decomposition product. The final dataset, now consisting of 2151 structures, was retrained by fine-tuning on the OMAT dataset with the multi-head fine-tuning approach, leading to a final model for each of the three interfaces of interest. Multi-head fine-tuning was enabled, and foundation_filter_element, foundation_model_readout, and foundation_model_elements were set to true. We used 100 000 samples of the pretrained head, a learning rate of 0.0001, and a weight decay of 1 × 10⁻⁹. The weights for energy, forces, and stresses were 1.0, 100.0, and 100.0 with no stage two employed. We trained the model for 30 epochs with a batch size of 8. Again, all other parameters were either kept at their default values or determined by the descriptors of the foundation model. The workflow is shown in Fig. 2.

Fig. 2 MACE model creation for each of the interface structures of interest. The active learning procedure is explained in detail in Fig. 1.

To validate the three MACE models, we calculated the root mean squared errors (RMSEs) of the energies and forces based on a test dataset consisting of 500 Structures taken from AIMD simulations of the respective interface at 400 K and 1000 K. The RMSEs are shown in Fig. S2 and were (100)/(100): 4.455 meV energies and 46.823 meV Å⁻¹ forces, (100)/(010): 3.126 meV energies and 35.864 meV Å⁻¹ forces, and (011)/(011): 3.122 meV energies and 44.669 meV Å⁻¹ forces. We also compared the average RDFs with AIMD results, shown in Fig. S3. The biggest discrepancy stems from the description of the S–P bond over the course of the trajectory, hence we also made sure the average S–P RDF fit the AIMD results to a satisfactory level (see Fig. S4).

2.3 Machine-learned MD

2.3.1 Interfaces. To create the larger supercells for the machine-learning MD (MLMD) simulations, we cut new surfaces from the bulk Li-argyrodite unit cell (52 atoms) with 6, 6, and 7 layers in the z-direction for the (100), (010), and (011) surfaces, respectively, and then created a 2 × 2 × 1 supercell. The same was done for the (100) and (011) Li-metal surfaces, both with 6 layers. One extra surface with the outer layer of the (100) surface removed was also created to ensure symmetry is conserved within the (100)/(100) interface supercell. We used INTERFACER to create the interface supercells from the existing surface structures, along with ASE's BFGS calculator with f_max = 0.04 meV Å⁻¹ to get an original optimisation of the atom positions within the interfaces. This was used as the starting structure for the MLMD simulations, which were run with LAMMPS.³⁶ We ran 10 ps NVT to equilibrate the system, using the Nose–Hoover thermostat with 1 fs timestep and a damping time parameter (t_damp) of 100 fs. After the equilibration, the ensemble was changed to NPT, again using Nose–Hoover with the same settings for the thermostat, and a target pressure of 1 bar and 1000 fs pressure damping factor (p_damp) for the barostat, with the anisotropic setting to control the 6 force components independently.

2.3.2 Bulk Li-argyrodite. Since the interface MLIP models were also trained on Li-argyrodite bulk structures at various temperatures, we decided to use one of the existing model for the MLMD simulations of the bulk rather than training a new one, after having ensured that they predict the bulk properties, such as the bulk RDF and coordination numbers, satisfactorily. To determine which of the three interface MACE models to use, we compared the RMSEs of the energies and forces of the three MLIPs on a test set made up of 400 Li-argyrodite bulk structures, with 100 structures from 400 K, 600 K, 800 K, and 1000 K, as shown in Fig. S5. The (100)/(010) model has the lowest errors for both energies and forces and was thus used for the MLMD runs. Setting up a 2 × 2 × 2 supercell of the bulk Li-argyrodite structure, we ran a 1 ns NPT simulation with LAMMPS after a 10 ps NVT equilibration. The settings for the NVT and NPT runs were the same as detailed for the interface structures above.

3 Results and discussions

3.1 Bulk structure

Li₆PS₅Cl crystallises in a cubic (F4̄3m) lattice, consisting of a face-centered arrangement of [PS₄]³⁻ polyhedra with P being on the 4b site and S on the tetrahedral 16e sites. The ‘free’ sulfur atoms that are not part of the [PS₄]³⁻ polyhedra and chlorine partially occupy the 4a and 4c sites. The Li ions are partially distributed across the 48h and 24g sites, forming cage-like polyhedra around the anion on the 4c site.

Earlier reports showed that having the S and Cl ions disordered across the two crystallographic sites is key for the superionic conductivity that can be observed in Li₆PS₅Cl compared to the ordered structure of Li₆PS₅I, where all the I is on the 4a sites and the free S on the 4c sites.^13,16,37,38 The exact extent of the disorder depends on the synthesis conditions but Cl distribution on 4a : 4c ranges between 38.5 : 61.5 and 46.2 : 53.8%.^38–42 To keep the computational cost low, we aimed to retain the original unit cell size with 4 formula units of Li₆PS₅Cl (52 atoms) and thus chose a 50 : 50 distribution across both sites. Both the ordered and disordered structures are shown in Fig. 3.

Fig. 3 (a) Ordered and (b) disordered structures of Li-argyrodite, with the 4a sites highlighted with circles and the 4c sites with squares. In the ordered structure, all Cl atoms are on the 4a sites, and all free S atoms on the 4c sites, while both species are distributed 50 : 50 across both sites in the disordered structure.

Geometry optimisation of the disordered structure led to a reduction of the F4̄3m symmetry to P1, indicating that the F4̄3m structure is metastable at 0 K. This is in line with what D'Amore et al.¹⁵ and Wang et al.⁴³ found previously, with the latter group showing that the cubic phase of Li₆PS₅Cl will be stable above 613.9 K. Our relaxed lattice parameters are in good agreement with the ones determined by D'Amore et al.¹⁵ and by Jeon et al.¹⁶ for the structure with a 50 : 50 distribution of S and Cl across the 4a and 4c sites (Table 1).

Table 1 Comparison of lattice constants of relaxed structures

Source	a (Å)	b (Å)	c (Å)	α (°)	β (°)	γ (°)
D'Amore^15	9.96	10.05	10.01	87.9	89.0	86.8
Jeon^16	9.71	9.83	9.71	89.0	89.3	89.9
This work	9.96	9.94	10.03	88.8	90.2	91.0

---

To compare the diffusion mechanism between the ordered and disordered structure, AIMD simulations were run in the NVT ensemble at 600 K for 40 ps for both models. Visualising the Li-ion trajectory during the simulation (see Fig. 4) shows the possible sites Li-ions can access and gives an indication of how often a given site is occupied. It can be seen that the anion-disordered structure has long-range diffusion, while in the ordered structure, the Li-ions only diffuse within the cages formed around the central S atoms.

Fig. 4 Li-ion trajectory lines of the (a) ordered and (b) disordered Li-argyrodite structures after 40 ps of AIMD simulation, with the starting positions of all the atoms in the system overlayed.

For a more detailed analysis of the diffusion behaviour of Li-ions in the disordered bulk structure and to obtain an accurate estimate of the diffusion coefficient and ionic conductivity, we trained an MLIP using MACE (see more details in Section 2.3.2). MD simulations in the NPT ensemble of a 2 × 2 × 2 supercell (containing 426 atoms) were then performed at 400 K, 600 K, 800 K, and 1000 K for 1 ns using the MLIP model.

The diffusion coefficient can be estimated from an MD trajectory using the Einstein relation:¹²

with D_s being the self-diffusion coefficient, N the total number of atoms of the diffusing species, and the mean squared displacement (MSD). Here we use kinisi,⁴⁴ which implemented an approximate Bayesian regression method to fit a linear model to the MSD data from the simulation. Since molecular dynamics (MD) data is neither statistically independent nor identically distributed, this gives a more accurate estimation of the diffusion coefficient and its related uncertainty than using ordinary least squares (OLS).⁴⁵ The MSD vs. time plots for each temperature are shown in Fig. S6.

The Nernst–Einstein equation can be used to determine the ionic conductivity from the self-diffusion coefficient when there is negligible correlation in ionic movement.¹²

where n is the ion density of Li, e the elementary electron charge, Z the valence of Li, k_B the Boltzmann constant, T the temperature, and D_s the self-diffusion coefficient. In the case of correlated ionic motion, as in our system, the charge diffusion coefficient D_σ, derived from the mean squared total displacement (MSTD) of the charge carriers' center-of-mass, needs to be used instead.⁴⁶ The plot of the MSTD vs. time at 400 K can be seen in Fig. S7.

While the room temperature diffusivity and conductivity are needed to compare with existing experimental results, the convergence of the trajectories is generally worse at lower temperatures due to slower particle motion. It is therefore preferred to run several trajectories at higher temperatures and extrapolate down to room temperature using the Arrhenius relation between D and 1/T, which also allows for the calculation of the activation energy (eV). This is more accurate than running a direct simulation at 300 K if the error of the diffusion coefficient obtained at each simulated temperature is kept low by ensuring the MSD is only calculated within the relevant region where the MSD vs. time (t) curve follows the linear Einstein relationship.^12,47

where D₀ is the maximum diffusivity at infinite temperature and E_a is the activation energy.

The diffusion coefficient of bulk Li-argyrodite at 300 K obtained from extrapolation is 5.79 × 10⁻⁷ ± 8.15 × 10⁻¹⁶ cm² s⁻¹, and the activation energy of the diffusion process E_a is 0.18 ± 5.14 × 10⁻⁶ eV (see Fig. 5). The same procedure was used to extrapolate the conductivity at 300 K, which was determined to be 1.39 × 10⁻¹ ± 1.73 × 10⁻⁴ S cm⁻¹. Due to the higher uncertainties when calculating the conductivity, the activation energy obtained from the Arrhenius plot of the conductivity (0.12 ± 1.91 × 10⁻⁵ eV) differs from the one obtained for the diffusion coefficient. The values in the literature for room temperature conductivity and the activation energy vary widely both for experimental and computational studies. Experimentally, the reported conductivities range from 10⁻² to 10⁻⁶ S cm⁻¹ and the activation energies lie between 0.11 to 0.45 eV. These differences stem from the processing methods and conditions, e.g. ball milling vs. solution-based processes as well as temperature effects, and the method used to measure the conductivity, e.g. impedance spectroscopy vs. ⁷Li NMR measurements.^{8,37–39,48–50} On the computational side, the room temperature conductivity varies between 6 × 10⁻⁶ and 0.89 S cm⁻¹ for the disordered system,^10,12,16 and the only reported activation energy for disordered systems we found is 0.25 eV.¹⁶ All these results were obtained from AIMD simulations of 100 ps to 300 ps using a simulation cell with 52 atoms. Our results are, to our knowledge, the only ones using an MLIP to go up to 1 ns simulations for a 2 × 2 × 2 supercell, which minimises finite size effects and maximises convergence, and sit well within the boundaries of the computational literature. Due to the absence of defects and grain boundaries, the determined value is however still higher than the experimentally determined conductivities.

Fig. 5 Arrhenius plot of σ vs. 1/T (orange) and D vs. 1/T (blue) of the bulk Li-argyrodite disordered structure, including the respective activation energies. Diffusion coefficients and conductivities were calculated at 400 K, 600 K, 800 K, and 1000 K.

3.2 Creation of interface models

3.2.1 Surface formation energies. To obtain a representative model of the Li-metal/Li-argyrodite interface, we first created various surface models and determined which ones are most likely to be stable. While the most stable surface will not necessarily form the most stable interface, it mimics the process where the electrode and the electrolyte are formed separately and subsequently pressed together under high pressures during cell manufacturing.^9,51,52 To limit computational efforts, we investigated all the surfaces identified with Miller planes up to 1. These include the planes that have a low lattice mismatch with Li-metal.¹⁷ Due to the complex crystal structure of Li-argyrodite, there are several possible surface terminations for each Miller plane, as shown in Fig. S8, which will affect the stability of the surface and the characteristics of the interface. To calculate the formation energy for each of them, it is important to symmetrise the surface models, as the resulting energy would otherwise be the average of two different surfaces. For each of the Miller planes, we created models of all possible symmetric surface terminations, leading to a total of 28 surface structures, and subsequently relaxed them with DFT. From the resulting DFT ground state energy, the surface formation energy (γ_surf), a measure by which to compare and rank the different surface models, can be calculated using the following equation:⁵³where A is the cross-sectional area of the surface (xy plane), E_surf the energy of the surface model as calculated by DFT, and E_bulk the energy of the bulk structure. Symmetrisation of the surface models leads to non-stoichiometric ratios. This is accounted for by n_i, the difference number between the surface and the bulk model for each chemical species i, and µ_i, the chemical potential of species i in Li-argyrodite.

The chemical potential of each species is determined by taking into account the decomposition products of Li-argyrodite, i.e. Li₃P, Li₂S, and LiCl.^18,54

E_Li6PS5Cl = 6µ_Li + µ_P + 5µ_S + µ_Cl (5)

E_Li3P = 3µ_Li + µ_P (6)

E_Li2S = 2µ_Li + µ_S (7)

E_LiCl = µ_Li + µ_Cl (8)

The γ_surf of all 28 surface structures are summarised in Table 2, Fig. S9 and S10.

Table 2 Surface formation energies (γ_surf) in eV nm⁻², chemical composition and charge of all surface models. The rows in bold indicate the three surfaces with the lowest surface formation energies

Miller index	Surface termination	γsurf (eV nm^−2)	Composition	Charge (e)
100	S(Li–Cl)	5.011	Li8PS7Cl1.5	−10
100	Li–Cl	2.203	Li8PS5Cl1.5	6
100	Li–P–S–Cl	1.006	Li10P2S8Cl1.5	10
100	S(Li–S)	5.691	Li8PS7.5Cl	−12
100	Li–S	3.033	Li8PS5.5Cl	4
100	Li–P–S–Cl	1.576	Li10P2S7.5Cl2	12
010	S(Li–S/Cl)	5.196	Li11P1.5S9.5Cl2	−10
010	Li–S/Cl	2.808	Li11P1.5S7.5Cl2	6
010	S(Li–P–Cl)	2.413	Li7P1.5S7Cl1.5	4
010	Li–P–Cl	1.225	Li7P1.5S5Cl1.5	12
010	S(Li–S/Cl)	5.632	Li11P1.5S10Cl1.5	−12
010	Li–S/Cl	2.795	Li11P1.5S8Cl1.5	4
010	S(Li–P–S)	3.208	Li7P1.5S7.5Cl	−6
010	Li–P–S	2.057	Li7P1.5S5.5Cl	10
001	S(Li–S/Cl)	5.448	Li8PS7.25Cl1.25	−11
001	Li–S/Cl	3.208	Li8PS5.25Cl1.25	5
001	Li–P–S–Cl	1.471	Li7P1.5S5.25Cl1.25	11
001	S(Li–P–S–Cl)	3.037	Li7P1.5S7.25Cl1.25	−5
110	Li–S	2.053	Li14P2S11Cl2	0
110	Li–P–S–Cl	1.861	Li16P3S14Cl3	1
101	Li–S	2.229	Li17P2.5S13.25Cl2.75	1
101	Li–P–S–Cl	2.174	Li13P2.5S11.25Cl2.75	1
101	Li–S	2.557	Li17P2.5S13.75Cl2.25	−1
101	Li–P–S–Cl	2.313	Li13P2.5S11.75Cl2.25	−1
011	Li–S	2.293	Li17P2.5S13.25Cl2.75	1
011	Li–P–S–Cl	1.168	Li13P2.5S11.25Cl2.75	1
011	Li–S	2.790	Li17P2.5S13.75Cl2.25	−1
011	Li–P–S–Cl	3.256	Li13P2.5S13.75Cl2.25	−1

---

Several trends can be observed from Table 2. The (100), (010), and (001) orientations are symmetrically equivalent for cubic systems, but the anion disorder in Li-argyrodite breaks that symmetry. As a result, there is a different number of symmetric terminations, even for these ‘related’ Miller indices, depending on the ordering of the free S and Cl atoms. Nevertheless, similar terminations give similar γ_surf for these three Miller planes. In contrast, no relation can be observed between surface terminations of the (110), (101), and (011) Miller planes. No (111) surfaces are included here, since there are no symmetric surface terminations for this Miller plane.

3.2.2 Interface formation energies. The three surfaces with the lowest γ_surf, i.e. the (100), (010), and (011) Li–P–S/Cl terminated models, were chosen to create interfaces with Li-metal. For the creation of the interface models, our in-house code INTERFACER was used. The lowest lattice mismatches were found for the (100)/(100), (100)/(010), and (011)/(011) Li-metal/Li-argyrodite interfaces, with all three having a mismatch of 0.37%. To ensure the interface models are also symmetric, the (100) Li-metal surface also had to be symmetrised by deleting one layer of Li-ion atom for the (100)/(100) interface. No changes were needed for the (100)/(010) interface since the resulting interfaces were identical but vertically translated, and the (011) surface of Li-metal is already symmetric.

For each model, the interface formation energy (γ_inter) was calculated according to Gao et al.⁵³

where A is the cross-sectional area of the interface (xy plane), E_total the DFT energy of the interface model, E_bulk1 and E_bulk2 the energies of the two bulk structures, n_i the difference number between the interface model and the bulk models for species i as explained above, and µ_i is the chemical potential of species i, as defined in eqn (5)–(8). By fitting both surfaces into one supercell, we induce an artificial strain onto the system that does not exist in a real physical system. Ideally, using adequately large models will minimise the lattice mismatch and hence the artificial strain, but this is not possible due to increased DFT costs. To approximate this strain, we relax the surface models with the a and b cell parameters fixed to the values of the relaxed interface model and subtract the energy of the freely relaxed surface, as introduced by Sen et al.⁵⁵

γ^min_inter = γ_inter − E^str (10)

where

E^str = E^str_slab1 + E^str_slab2 (11)

andwith E^strained_slab and E^relaxed_slab being the energy of the strained and relaxed slab, respectively, and A being the surface area. The resulting strain values are summarised in Table 3.

Table 3 E_str in eV nm⁻² for the different Li-argyrodite and Li-metal surfaces

Miller index	Bulk	Estr (eV nm^−2)
(100)	Li-argyrodite	−0.074
(010)	Li-argyrodite	0.802
(011)	Li-argyrodite	0.618
(100)	Li-metal	0.011
(100) sym	Li-metal	0.229
(011)	Li-metal	0.417

---

3.2.3 Charge-balancing. When symmetrising the surfaces, the resulting surface slabs are no longer stoichiometric, inducing an artificial charge. The effect of this surface charge could be compensated by having a large enough slab to imitate the infinite reservoir of electrons from the bulk that a real material surface would have available. However, prohibitively many atoms would be required for such a model. We thus take an approach similar to Canepa et al.,⁵⁶ where we delete the appropriate number of atoms to balance the charge, opting to always delete the minimum number of atoms to affect the structure as little as possible. For example, the chosen (100) surface with Li–P–S–Cl termination has the chemical formula of Li₄₀P₈S₃₂Cl₆, and thus a charge of +10, with the formal charges being Li: +1, P: +5, S: −2, and Cl: −1. Using configuration enumeration, we determine the two phosphorous atoms which are energetically the most favourable to remove from the surface (see Fig. S11). The same was done for the other two low energy surfaces, and the determined atoms were also removed from the three interface models.

The charge-balanced interface structures before and after relaxation are shown in Fig. 6 with their respective formation energies. The final models contain 99 Li-metal/84 Li-argyrodite atoms for the (100)/(100) interface, 108 Li-metal/56 Li-argyodite atoms for the (100)/(010) interface, and 108 Li-metal/117 Li-argyrodite atoms for the (011)/(011) interface. The changes in γ_surf and γ_inter between the charged and the charge-balanced models are summarised in Table 4, and the interface structures prior to charge-balancing are shown in Fig. S12. While the changes in γ_surf correlate with the amount of charge, the same is not true for the γ_inter. The effect of the charge-balancing on the interface depends on how the removal of the atom(s) affect the interaction with Li-metal as well, and while γ_inter of the (100)/(100) and (100)/(010) interfaces become less negative, it becomes more negative for the (011)/(011) interface, where only a single Li atom is removed.

Fig. 6 (a) (100)/(100), (b) (100)/(010), and (c) (011)/(011) Li-metal/Li-argyrodite interfaces before (top) and after (bottom) geometry optimisation with DFT after charge-balancing. The Miller plane for Li-metal was chosen so that the interface would have minimal lattice mismatch. The interface formation energies (γ_inter) are given for each structure in eV nm².

Table 4 Surface formation energy (γ_surf) and interface formation energy (γ_inter) of charged and neutral surface and interface models in eV nm⁻²

γsurf (eV nm^−2)	(100)	(010)	(011)
Charged	1.006	1.225	1.168
Neutral	3.039	2.834	1.261

---

γinter (eV nm^−2)	(100)/(100)	(100)/(010)	(011)/(011)
Charged	−12.903	−11.835	−8.410
Neutral	−7.807	−8.584	−12.947

---

3.3 Analysis of interface models

3.3.1 Chemical stability of the Li-metal/Li-argyrodite interface. While no formation energies have been reported for Li-metal/Li-argyrodite interfaces, both Camacho-Foreroa et al.⁵⁷ and Lepley et al.⁵⁸ have found negative formation energies for Li₃PS₄, a decomposition product of Li-argyrodite. They argued that this is due to the adhesive forces between the materials being stronger than the cohesive forces within the individual materials, thus showing that reactions are spontaneously happening at the interface, indicating an unstable interface.

Further analysis of the optimised structures also shows that Li-argyrodite is not stable against Li-metal, in line with previous reports.^9,17,18,59 From visual inspection, it can be seen that an amorphous layer forms at the interface during optimisation. The RDFs of the P–S and the Li–P bonds of the (100)/(100) interface show a decrease in the P–S peak intensity while new P–Li bonds at various new bond lengths form (see Fig. 7a). This indicates that the PS₄ polyhedra break down and P forms Li_xP species with Li, as also observed by Cheng et al.¹⁸ for the Li (100)/Li₆PS₅Cl (100) interface.

Fig. 7 (a) Radial distribution functions of the P–S (top) and P–Li (bottom) bond lengths of the (100)/(100) interface before and after the geometry optimisation. (b) DDEC6 charge partitioning of the optimised (100)/(100) interface.

The same decomposition reactions can be observed when analysing the partial charges, obtained through DDEC6 charge partitioning²⁸ (see Fig. 7b). When the PS₄ polyhedra break down and P reacts with Li-metal from the interface, P changes its atomic charge from ∼0.8 to ∼−2. These reactions lead to a reduction in the total number of metallic Li, and the S of the broken-down polyhedra also change their atomic charge from ∼−0.8 to ∼−1.6, similar to the free S atoms. There appears to be a correlation between the number of broken down PS₄ polyhedra and the magnitude of the negative interface formation energy. These spontaneous reactions can also be seen in the other two interfaces (see Fig. S13) during DFT geometry optimisations at 0 K, indicating a kinetically barrier-less process. To probe performance of Li-argyrodite as solid electrolyte when coupled with Li-metal, it is important to investigate the dynamic properties of the SEI layer formed.

For a detailed analysis of the stability of the interface models and to assess their suitability for use in ASSBs, we built an MLIP model for each interface under consideration. This enabled us to increase the size of the interface cells to 2220, 2352, and 2596 atoms for the (100)/(100), (100)/(010), and (011)/(011) interfaces, respectively, and run MD simulations of 5 ns at 300 K and 400 K, enabling more realistic insights into the processes happening at the interfaces.

The formation and growth of the SEI layer can be seen clearly when visualising the 5 ns MD trajectories at 400 K (see Fig. 8 and S14–S19). All three interfaces decompose and form an originally amorphous layer, which starts to crystallise within the first 1–2 ns. The crystallisation happens at one of the two Li-metal/Li-argyrodite interfaces present within each model and moves towards the other interface until the whole SEI is a crystalline layer. The (011)/(011) interface only fully crystallises towards the end of the 5 ns. To ensure that the results are converged and no further changes happen, we ran the simulation for another 5 ns, resulting in a total of 10 ns. Due to the slower particle movement at 300 K, the decomposition proceeds at least 5 times slower, to the point that the decomposition is not completed at 5 ns in the (100)/(010) interface. The onset of a crystalline region can only be observed in the (100)/(010) interface (see Fig. S18), and no crystalline phase forms within 5 ns in the other two interfaces. The exact times of decomposition and crystallisation are given in Table 5.

Fig. 8 Snapshots of the 5 ns MD trajectory of the (100)/(100) interface at 400 K, taken every nanosecond. The vertical thick dashed lines indicate the position of the interface, determined by the non-Li atom that diffused furthest into the Li-metal. The thin lines show the crystalline regime which grows throughout the simulation. The red dotted lines at t = 0 ns indicate the position of the membrane used to calculate the Flux in the regions (1) interface, (2) electrode, and (3) electrolyte. The membrane position at the interface is placed at the actual interface position at time 0.

Table 5 End time of the decomposition (D), as well as start and end time of crystallisation (C) for the (100)/(100), (100)/(010), and (011)/(011) interface at 300 K and 400 K

Interface	T (K)	End D (ns)	Start C (ns)	End C (ns)
(100)/(100)	300	1.93	—	—
	400	0.39	0.24	2.43
(100)/(100) sym	400	0.35	4.73	>5
(100)/(010)	300	>5	1.76	>5
	400	0.57	0.87	1.9
(011)/(011)	300	2.11	—	—
	400	0.37	0.43	4.80

---

The changes in the coordination environment during the MD simulations were analysed by tracking the Li–P, Li–S, Li–Cl, and P–S bonds (Fig. 9 and S20). Within the first 500 ps the Li–S and Li–P coordination numbers increase to around 8 and 10, respectively, while the S–P coordination number goes down to 0. This shows that all the PS₄ polyhedra in the system break down and Li₂S (CN = 8) and Li₃P (CN = 11) as well as other Li_xP species form, making up the amorphous layer together with LiCl. The coordination number of Li₃P starts decreasing after reaching its maximum at around 500 ps, while the LiCl coordination number increases, showing that P and Cl are in a different coordination environment in the crystalline phase compared to the amorphous phase. Closer analysis of the coordination environment reveals that the phase crystallises in the antifluorite structure. By classifying all Li atoms that have a coordination number of 3 or more to S, P, or Cl as part of the SEI, we can determine the composition in the (100)/(100) system as Li₂S_0.72P_0.14Cl_0.14. The crystalline phase can thus be described as Li₂S, with S atoms being replaced by P and Cl. Using the same method, the exact composition varies slightly for the (100)/(010) and (011)/(011) interfaces, with Li₂S_0.71P_0.14Cl_0.15 and Li₂S_0.70P_0.15Cl_0.15. These small variations can be due to defects within the crystal structure and the simple way to classify Li-metal ions. Overall it can be said that the crystalline phase that forms is the same in all three systems, as shown in Fig. 10. While it appears that the (011)/(011) interface has not fully crystallised at 5 ns, the system actually forms a polycrystalline SEI, with the final part of SEI on the left of the Li-metal slab crystallising in a different orientation (see Fig. S18). Analysis of the composition of this phase and the coordination numbers shows that the second phase also crystallises in an antifluorite arrangement, as is shown in Fig. S21.

Fig. 9 Coordination number analysis of S–P, S–Li, P–Li, and Cl–Li bonds for the (100)/(100) interface using MLMD trajectories at 300 K and 400 K.

Fig. 10 Crystal structure of the (a) (100)/(100), (b) (100)/(010), and (c) (011)/(011) interfaces at 5 ns. Note that due to the change in perspective to capture the antifluorite structure, the edges of (011)/(011) SEI regions overlap, which does however not represent a full penetration of the Li-metal region by the SEI.

This matches what Chaney et al.²⁰ and Ding et al.²¹ reported. However, in contrast to our results, both groups observed the growth of the crystalline phase to simultaneously start and proceed from both interfaces. One possible reason for this difference could be that while we initially generated perfectly symmetric surface cuts, we did not preserve symmetry when charge-balancing, but instead removed the energetically most favourable atoms. To test this hypothesis, we created a fully symmetric supercell for the (100)/(100) interface, where it was possible to remove the atoms in a way that keeps the symmetry of the system intact, and ran 5 ns of MD at 400 K using the existing MLIP model for the (100)/(100) interface. The symmetric system only starts to crystallise towards the end of the 5 ns simulation, with the crystallisation starting at the centre of the amorphous SEI region. This confirms that the exact charge environment of the interface does have an effect on the rates of decomposition and crystallisation, and the charge imbalance leads to enhanced nucleation.

Furthermore, we can see that the surface termination also has an effect on the stability of the Li-metal. Li only stays in its BCC arrangement in the (100)/(100) interface at both 300 K and 400 K, while it melts during the 400 K (100)/(010) and the 300 K and 400 K (011)/(011) MD simulations, despite the melting points of Li being 453.65 K. The melting can also be observed in 20 ps AIMD simulations of the smaller interface models, showing that it is not due to a failure of the MLIP model. The unexpected melting behaviour is likely a result of local pressure build-up in the NVT ensemble, a reduction of the melting point when the anionic species penetrate into the Li-metal layer and a new amorphous layer is formed, and/or finite-size effects stemming from the size limitations of ab initio methods. Similar behaviour can be observed in AIMD simulations of interfaces between Li-metal and thiophosphate SSEs in previously published work.^17,19,57 The melting of Li-metal does, however, not appear to affect the onset or speed of the crystallisation.

To investigate the rate of decomposition and the phase transition from amorphous to crystalline, we calculated the thickness of the SEI layer and the degree of crystallisation. For the SEI layer thickness we are using the same approach as Ding et al.,²¹ determining the number of intact PS₄ polyhedra to classify the region of Li-argyrodite, and calculating the distance between the outer edges of this region and the most outer non-Li atom (S, P or Cl) at the interface to Li-metal. The degree of crystallisation (in %) was determined by counting all the S atoms with a coordination number of 8 that would be expected in the antifluorite structure, and dividing this by the total number of S atoms in the system. Looking at the (100)/(100) system at 400 K, the SEI thickness increases rapidly to ∼64 Å at 500 ps and then continues increasing to its maximum of ∼66 Å at 1 ns at a slower rate, before it decreases slightly and reaches a fixed thickness of ∼64 Å (Fig. 11a). The decrease can be associated with the difference in volume between the amorphous and the crystalline phase, and is more pronounced for the other two interfaces (see Fig. S22). At 300 K, the SEI thickness of the (100)/(100) interface increases during the first 500 ps up to ∼62 Å, where it stabilises. This confirms that the decrease that can be seen in the 400 K trajectories is due to the crystallisation, which does not start in 5 ns for the (100)/(100) interface at 300 K.

Fig. 11 (a) Evolution of the SEI thickness (in Å) of the (100)/(100) interface, calculated using the distance between the most outer atom of Li-argyrodite with intact PS₄ polyhedra, and the atom that diffused furthest into Li-metal. (b) Evolution of the degree of crystallisation of the SEI in % in the (100)/(100) interface, calculated by determining the number of S with 8 bonds to Li, as in Li₂S, and dividing it by the total number of S atoms in the system.

At 400 K, the crystallisation has an equally rapid increase within the first 500 ps up to ∼55%, before it increases more slowly to its maximum value just below 80% during the phase transition from amorphous to crystalline. The same trend can be seen for the (100)/(010) and (011)/(011) 400 K simulations (Fig. S23). The crystallisation does not reach 100% due to grain boundaries and defects in the system, as can be seen in Fig. 10.

From the 300 K trajectory, it can be seen that the degree of crystallisation stays around 50–55% in the amorphous layer, so only values above 55% indicate the presence of the crystalline region. This is because Li₂S is present as individual, disconnected species in the amorphous arrangement, while long-range order only forms during the crystallisation. The way the crystallisation is calculated does, however, only calculate how many Li–S bonds exists and not explicitly capture how many Li₂S particles are connected.

To verify that the dynamics of Li-argyrodite decomposition into an amorphous phase and subsequent crystallisation are reproducible and not a result of the specific atoms chosen for deletion to balance the charges in the system, we ran 3 ns simulations of the (100)/(100) interface of two systems with different P atoms removed from the surface at 400 K. 3 ns was chosen as the original system is fully crystallised at this point and no further changes can be observed afterwards. The changes in coordination numbers are shown in Fig. S24, and the evolution of the SEI thickness and the crystallisation fraction are compared to the original model in Fig. S25. While there are differences in how fast the crystallisation proceeds, the point of full decomposition and onset of crystallisation are sufficiently similar (see Table S1). Furthermore, we confirmed that the SEI of the other two models also crystallise in the antifluorite phase, as shown in Fig. S26. This gives us confidence that the results are not artificially biased and the simulations are general and reproducible.

3.3.2 Diffusion properties of the Li-metal/Li-argyrodite interfaces. To better understand how the formation of the SEI affects the ion transport, and hence the performance of the battery, we studied the Li-ion diffusion across the interface and within the SEI layer. Noting that Li atoms have distinct diffusion behaviour in different regions of the battery system, e.g. bulk Li-metal, SEI layer and bulk Li-argyrodite, one should analyse them separately rather than using an average Li diffusion for the whole model. Obtaining the diffusion coefficient via the MSD to this end, as was done for bulk Li-argyrodite, is not straightforward, since many individual Li-atoms cross into a different region during the simulations. One needs to constantly track the type of each Li-ion and keep an account of their contribution to each MSD variable defined, making this a computationally arduous task.

A possible solution would be to calculate MSDs over short enough time periods that no diffusion into other regions is observed, but this would prevent the MSD convergence and thus result in diffusion coefficients that are not statistically meaningful. Instead, we calculate the bi-directional flux of Li atoms across a 3 Å-thick membrane, similar to Sen et al.,⁵⁵ which can be placed anywhere within the simulation cell. We investigated the interface region at the centre of the cell, determined by the non-Li atom that has diffused into Li-metal the furthest during the equilibration run. This corresponds to the position of the dashed red line at 0 ns in Fig. 8 and is labelled with (1). The relative membrane position within the cell has to be kept constant for the calculation of the flux to be accurate (their positions are thus defined using fractional coordinates). However, as can be seen from the other snapshots in Fig. 8, the decomposition and consequent reduction of Li-metal leads to the shrinkage of the bulk Li-metal region and a shift of the interface towards the left, and the membrane will thus eventually be located in the SEI layer.

Fig. 12a shows that the flux of the (100)/(100) interface at 400 K from left to right, i.e. from Li-metal towards Li-argyrodite, is larger than the other way around, with the only driving force being the differences in chemical potentials between the regions. The flux from left to right can be separated into three different regions: (1) quickly rising flux for the first 500 ps corresponding to the decomposition, (2) increasing flux with a less steep slope up to 2.2 ns, and (3) no flux with slope levelling off for the remainder of the 5 ns simulation. The start of region (2) and (3) coincide with the onset and completion of the crystallisation (see Table 5). The flux from right to left in contrast has only two distinguishable regions, one with constant increase in flux and one without any increase in flux, with region 2 also starting at around 2.2 ns. For the (100)/(010) and the (011)/(011) interfaces, shown in Fig. S27a and S28a, the different regions are less distinguishable and the flux grows exponentially to start with. The (011)/(011) interface does not reach a plateau, which is most likely due to the membrane position being closer to the actual interface than for the other two systems, leading to the movement of Li-ions close to the interface being captured.

Fig. 12 Flux from left to right (L → R) and right to left (R → L) at the centre of the interface, electrode and electrolyte of the (100)/(100) interface at 300 K and 400 K, calculated by counting the number of atoms moving across a 3 Å-thick membrane.

Moving the membrane towards the centre of the Li-metal, indicated by the dotted red line at position (2) in Fig. 8, allows us to explore the flux in the electrode. For the (100)/(100) interface (see Fig. 12b), there is minimal (400 K) to no (300 K) flux, showing that the Li-metal stays within the BCC arrangement, where there is no long-range movement of Li-atoms. As mentioned above, the other two interfaces experience melting of the Li-metal, which can be observed in the flux and is shown in Fig. S27b and S28b. At 400 K the (100)/(010) interface only stays stable for the first 500 ps before melting, shown by the steep increase in the flux, while at 300 K the Li-metal is stable and only a minimal flux of Li-atoms is observed. For the (011)/(011) interface, melting is observed at both temperatures. At 400 K the system melts instantly, with a decrease in the slope of the flux after 1.7 ns, while at 300 K the Li-metal only melts after 1.5 ns.

The last region we explored is the centre of the bulk electrolyte, highlighted by the dotted red line at position (3) in Fig. 8. For the (100)/(100) interface at 400 K, shown in Fig. 12c, the flux increases until 1.2 ns, where it levels off. The flux from left to right again has a slightly larger slope than from right to left. This is the same for the (100)/(010) system, as can be seen in Fig. S27c. The flux of the (011)/(011) interface (see Fig. S28c) follows the same trend to start with, but does not converge after the crystallisation of the SEI is completed and instead continues to increase with a small slope, especially from left to right. The difference in flux from left to right and right to left indicates a driving force towards the right, likely due to the difference in chemical potential at each interface that was introduced by the non-symmetric deletion of P and Li atoms for charge-balancing.

At 300 K, the flux from right to left and left to right are nearly equal, and roughly follow a dependence, although with some small changes in the slope. Similar behaviour can be observed for the (100)/(010) and (011)/(011) interfaces. However, the flux from left to right and right to left are not equal in the latter two systems. For the (100)/(010) interface, the flux from right to left is larger, while in the (011)/(011) system, the flux from left to right dominates. The exact values of the flux at each membrane position in both directions at 300 K and 400 K for each of the three interfaces are summarised in Table S2.

All three interface systems shrink initially and stabilise after a few ns. To obtain the flux at the interface after the volume stabilised, and of the crystalline SEI, we calculated the bi-directional flux for one ns starting from 4 ns of the trajectories of all three interfaces. To obtain the flux of the amorphous SEI, we placed a new membrane at one quarter of the c- of the electrolyte for the (100)/(100) and (100)/(010) interfaces, where the SEI grows from left to right, and three quarters for the (011)/(011) interface, where the SEI grows in the other direction. We calculated the flux at these new membranes for a 1 ns time interval where the membrane was continuously in a region of amorphous SEI. The time intervals chosen were 1–2 ns for the (100)/(100) and (011)/(011) interfaces, and 0.5–1.5 ns for the (100)/(010) interface. The membrane positions are defined with respect to the cell at the start point of that interval. The results are summarised in Table 6.

Table 6 Average flux at each interface and the center of the amorphous and crystalline SEI for one ns of the trajectory of the (100)/(100), (100)/(010), and (011)/(011) interfaces at 400 K

Interface	Position	Flux L → R (atoms/nm^2 ns)	Flux R → L (atoms/nm^2 ns)
(100)/(100)	Interface	65.25	64.12
	Amorphous SEI	12.94	12.26
	Crystalline SEI	0.23	0.23
(100)/(010)	Interface	79.57	78.92
	Amorphous SEI	10.41	18.22
	Crystalline SEI	0.0	0.0
(011)/(011)	Interface	67.88	66.96
	Amorphous SEI	6.82	13.79
	Crystalline SEI	0.93	0.46

---

The crystalline SEI has the highest flux in the (011)/(011) interface, explaining why the flux at the centre of the SEI does not fully level off even after the crystallisation has happened. At 400 K, the average flux of bulk Li-argyrodite from left to right and left to right is 35.142 atoms/nm² ns and 37.032 atoms/nm² ns, respectively. Fig. 13 highlights the difference in the total flux of the bulk system, the amorphous SEI, and the crystalline SEI for the (100)/(100) interface. The same trend is seen for the other two interfaces (see Fig. S29). This shows that the flux at the actual interface is not the bottleneck, but that instead the low flux of the crystallised SEI is the inhibiting factor.

Fig. 13 Total flux in atoms per nm² per ns of bulk Li-argyrodite along with the amorphous and crystalline SEI for the (100)/(100) interface, tracked for 1 ns.

In contrast to this, we have previously shown that Li₃P–Li₂S solid solutions, which crystalise in the Li₂S anti-fluorite phase in a fully-reduced form, exhibit conductivities comparable to other SSE materials, while being chemically stable against Li-metal (no degradation behaviour).^55,60 The difference is most likely due to the presence of chlorine atoms in the Li-argyrodite SEI layer, less able to coordinate Li atoms due to its less favourable electrostatic interactions, and the high-energy interstitial sites needed for percolated migration can thus not be activated as easily.

4 Conclusions

We created three representative models of the Li-metal/Li-argyrodite interface by identifying three Li-argyrodite surfaces with the lowest formation energies, yielding interfaces with minimal lattice mismatch with Li-metal, namely Li-metal/Li-argyrodite (100)/(100), (100)/(010), and (011)/(011) interfaces. By training custom-made MACE-based MLIP models, we studied the dynamic processes at these interfaces through 5 ns MD simulations of interface systems with over 2000 atoms. The MLMD simulations at 300 K and 400 K revealed the initial growth of the SEI layer, which forms due to the chemical instability of Li-argyrodite against Li-metal. Initially, an amorphous phase containing the decomposition products of Li₃P, Li₂S and LiCl forms for all three interface models. At 400 K, this layer then crystallises into an antifluorite structure, with a composition of Li₂S_1−x−yP_xCl_y; {x, y = 0.14–0.15}. The thickness of the SEI decreases slightly during the crystallisation due to the more close-packed atomic arrangement in the crystalline phase. In all three interfaces, we observe complete decomposition of the Li-argyrodite within 5 ns at 300 K and 400 K, and the whole SEI layer crystallises within the 5 ns at 400 K in all three systems.

Calculating the Li-ion flux in various regions of the interface systems, namely the interface, bulk electrode (Li-metal), and bulk electrolyte (Li-argyrodite) regions, we showed that the flux decreases during the crystallisation before it completely levels off. The flux within the crystalline SEI phase drops to near zero, indicating that while ion transport is still possible in the amorphous layer, the Li-ions can no longer move freely in the crystalline layer. This highlights the high interfacial resistance stemming from the crystallisation of the amorphous SEI layer, which completely blocks the Li-ion diffusion. While the SEI layer formed is very similar for all three interfaces simulated in terms of crystal structure and chemical composition, there are differences in the rates of decomposition and crystallisation. The exact chemical composition and charge environment at the interface appear to affect the onset of nucleation of the crystalline phase. These atomic-level insights could potentially be used to control the SEI formation in sulphide-based SSEs and others, in designing approaches to prevent crystallisation of SEI layers.

Due to the limitations of the atomistic scale simulations, and currently poor description of long-range effects and electrostatics of machine-learning models,⁶¹ this work only considers a half cell at first contact, with the difference in chemical potentials being the only driving force for Li diffusion. Many chemo-mechanical phenomena happening in the full systems during Li stripping and plating, such as the formation of voids at the Li-metal/Li-argyrodite interface, which can lead to dendrite growth at high enough current densities,⁷ can therefore not be captured by our models. Lots of work is being put into developing MLIPs that can accurately describe these effects, hopefully enabling the simulation of the full ASSB under cycling conditions in future.^62,63

Author contributions

Chantal Baer: conceptualization, data curation, formal analysis, investigation, methodology, validation, visualization, writing – original draft. Roman Shantsila: software, writing – review & editing. Lukasz Figiel: conceptualization, funding acquisition, project administration, resources, supervision, writing – review & editing. Bora Karasulu: conceptualization, funding acquisition, methodology, project administration, resources, software, supervision, writing – review & editing.

Conflicts of interest

There are no conflicts to declare.

Data availability

All relevant bulk, surface, and interface structures used for this work are available on GitHub: https://github.com/bkarasulu/Li-argyrodite-interface. Supplementary information: plots and tables for model validation and extra results. See DOI: https://doi.org/10.1039/d6ta00922k.

Acknowledgements

C. B. and R. S. acknowledge the UK Engineering and Physical Sciences Research Council (EPSRC) funded CDT for Modelling of Heterogeneous Systems (HetSys CDT) [EP/S022848/1] for the PhD studentships. B. K. acknowledges support from the 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. We also acknowledge computational resources from ARCHER2 UK National Computing Service which was granted via ‘Access to HPC’ call project (e897).

References

1. T. Famprikis, P. Canepa, J. A. Dawson, M. S. Islam and C. Masquelier, Nat. Mater., 2019, 18, 1278–1291.
2. J. Janek and W. G. Zeier, Nat. Energy, 2016, 1, 16141.
3. S. A. Pervez, M. A. Cambaz, V. Thangadurai and M. Fichtner, ACS Appl. Mater. Interfaces, 2019, 11, 22029–22050.
4. D. Lin, Y. Liu and Y. Cui, Nat. Nanotechnol., 2017, 12, 194–206.
5. L. Sultanova and L. Figiel, Comput. Mater. Sci., 2021, 186, 109990.
6. S. Bazazzadeh, M. Pasta and L. Figiel, Int. J. Mech. Sci., 2024, 264, 108808.
7. J. Kasemchainan, S. Zekoll, D. S. Jolly, Z. Ning, G. O. Hartley, J. Marrow and P. G. Bruce, Nat. Mater., 2019, 18, 1105–1111.
8. H. J. Deiseroth, S. T. Kong, H. Eckert, J. Vannahme, C. Reiner, T. Zaiß and M. Schlosser, Angew. Chem., Int. Ed., 2008, 47, 755–758.
9. C. Yu, L. van Eijck, S. Ganapathy and M. Wagemaker, Electrochim. Acta, 2016, 215, 93–99.
10. N. J. J. de Klerk, I. Rosłoń and M. Wagemaker, Chem. Mater., 2016, 28, 7955–7963.
11. Z. Deng, Z. Zhu, I. H. Chu and S. P. Ong, Chem. Mater., 2017, 29, 281–288.
12. A. Baktash, J. C. Reid, T. Roman and D. J. Searles, npj Comput. Mater., 2020, 6, 162.
13. B. J. Morgan, Chem. Mater., 2021, 33, 2004–2018.
14. N. T. Han, W. Bang-Li, K. I. Lin, V. K. Dien and M. F. Lin, RSC Adv., 2022, 12, 32674–32683.
15. M. D'Amore, L. E. Daga, R. Rocca, M. F. Sgroi, N. L. Marana, S. M. Casassa, L. Maschio and A. M. Ferrari, Phys. Chem. Chem. Phys., 2022, 24, 22978–22986.
16. T. Jeon, G. H. Cha and S. C. Jung, J. Mater. Chem. A, 2024, 12, 993–1002.
17. A. Golov and J. Carrasco, ACS Appl. Mater. Interfaces, 2021, 13, 43734–43745.
18. T. Cheng, B. V. Merinov, S. Morozov and W. A. Goddard, ACS Energy Lett., 2017, 2, 1454–1459.
19. J. Wang, A. A. Panchal, G. S. Gautam and P. Canepa, J. Mater. Chem. A, 2022, 10, 19732–19742.
20. G. Chaney, A. Golov, A. van Roekeghem, J. Carrasco and N. Mingo, ACS Appl. Mater. Interfaces, 2024, 16, 24624–24630.
21. J. Ding, L. Zichi, M. Carli, M. Wang, A. Musaelian, Y. Xie and B. Kozinsky, Coupled reaction and diffusion governing interface evolution in solid-state batteries, arXiv, 2025, preprint, arXiv:2506.10944,  DOI:10.48550/arxiv.2506.10944.
22. A. Jain, S. P. Ong, G. Hautier, W. Chen, W. D. Richards, S. Dacek, S. Cholia, D. Gunter, D. Skinner, G. Ceder and K. A. Persson, APL Mater., 2013, 1, 011002.
23. G. Kresse and J. Hafner, Phys. Rev. B, 1993, 47, 558–561.
24. G. Kresse and J. Furthmüller, Phys. Rev. B, 1996, 54, 11169–11186.
25. G. Kresse and J. Furthmüller, Comput. Mater. Sci., 1996, 6, 15–50.
26. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1997, 78, 1396.
27. G. Kresse and D. Joubert, Phys. Rev. B, 1999, 59, 1758–1775.
28. T. A. Manz and N. G. Limas, RSC Adv., 2016, 6, 47771–47801.
29. I. Batatia, D. P. Kovacs, G. N. C. Simm, C. Ortner and G. Csanyi, Adv. Neural Inf. Process. Syst., 2022.
30. I. Batatia, S. Batzner, D. P. Kovács, A. Musaelian, G. N. C. Simm, R. Drautz, C. Ortner, B. Kozinsky and G. Csányi, The Design Space of E(3)-Equivariant Atom-Centered Interatomic Potentials, 2022.
31. R. Shantsila, C. M. I. B. Baer, A. Bartók-Pártay and B. Karasulu, ChemRxiv, 2026, preprint,  DOI:10.26434/chemrxiv.15001544/v1.
32. I. Batatia, et al., J. Chem. Phys., 2025, 163, 184110.
33. L. Barroso-Luque, M. Shuaibi, X. Fu, B. M. Wood, M. Dzamba, M. Gao, A. Rizvi, C. L. Zitnick and Z. W. Ulissi, Open Materials 2024 (OMat24) Inorganic Materials Dataset and Models, arXiv, 2024, preprint, arXiv:2410.12771,  DOI:10.48550/arxiv.2410.12771.
34. A. Hjorth Larsen, et al., J. Phys.: Condens. Matter, 2017, 29, 273002.
35. M. Radova, W. G. Stark, C. S. Allen, R. J. Maurer and A. P. Bartók, npj Comput. Mater., 2025, 11, 237.
36. 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.
37. I. Hanghofer, M. Brinek, S. L. Eisbacher, B. Bitschnau, M. Volck, V. Hennige, I. Hanzu, D. Rettenwander and H. M. Wilkening, Phys. Chem. Chem. Phys., 2019, 21, 8489–8507.
38. P. R. Rayavarapu, N. Sharma, V. K. Peterson and S. Adams, J. Solid State Electrochem., 2012, 1807–1813.
39. M. A. Kraft, S. P. Culver, M. Calderon, F. Böcher, T. Krauskopf, A. Senyshyn, C. Dietrich, A. Zevalkink, J. Janek and W. G. Zeier, J. Am. Chem. Soc., 2017, 139, 10909–10918.
40. R. Schlenker, A. L. Hansen, A. Senyshyn, T. Zinkevich, M. Knapp, T. Hupfer, H. Ehrenberg and S. Indris, Chem. Mater., 2020, 32, 8420–8430.
41. N. Minafra, M. A. Kraft, T. Bernges, C. Li, R. Schlem, B. J. Morgan and W. G. Zeier, Inorg. Chem., 2020, 59, 11009–11019.
42. C. Liu, B. Chen, T. Zhang, J. Zhang, R. Wang, J. Zheng, Q. Mao and X. Liu, Angew. Chem., Int. Ed., 2023, 62, e202302655.
43. Z. Wang and G. Shao, J. Mater. Chem. A, 2017, 5, 21846–21857.
44. A. R. McCluskey, A. G. Squires, J. Dunn, S. W. Coles and B. J. Morgan, J. Open Source Softw., 2024, 9, 5984.
45. A. R. McCluskey, S. W. Coles and B. J. Morgan, J. Chem. Theory Comput., 2025, 21, 79–87.
46. A. Marcolongo and N. Marzari, Phys. Rev. Mater., 2017, 1, 025402.
47. X. He, Y. Zhu, A. Epstein and Y. Mo, npj Comput. Mater., 2018, 4, 18.
48. N. C. Rosero-Navarro, A. Miura and K. Tadanaga, J. Sol-Gel Sci. Technol., 2019, 89, 303–309.
49. D. H. Kim, D. Y. Oh, K. H. Park, Y. E. Choi, Y. J. Nam, H. A. Lee, S. M. Lee and Y. S. Jung, Nano Lett., 2017, 17, 3013–3020.
50. H. J. Deiseroth, J. Maier, K. Weichert, V. Nickel, S. T. Kong and C. Reiner, Z. Anorg. Allg. Chem., 2011, 637, 1287–1294.
51. Y. G. Lee, S. Fujiki, C. Jung, N. Suzuki, N. Yashiro, R. Omoda, D. S. Ko, T. Shiratsuchi, T. Sugimoto, S. Ryu, J. H. Ku, T. Watanabe, Y. Park, Y. Aihara, D. Im and I. T. Han, Nat. Energy, 2020, 5, 299–308.
52. Y. J. Choi, Y. J. Hwang, S. I. Kim and T. Kim, ACS Appl. Energy Mater., 2024, 7, 4421–4428.
53. B. Gao, R. Jalem and Y. Tateyama, ACS Appl. Mater. Interfaces, 2020, 12, 16350–16358.
54. C. Yu, F. Zhao, J. Luo, L. Zhang and X. Sun, Nano Energy, 2021, 83, 105858.
55. H. S. Sen and B. Karasulu, J. Mater. Chem. A, 2025, 13, 19878–19895.
56. P. Canepa, J. A. Dawson, G. S. Gautam, J. M. Statham, S. C. Parker and M. S. Islam, Chem. Mater., 2018, 30, 3019–3027.
57. L. E. Camacho-Foreroa and P. B. Balbuenaa, J. Power Sources, 2018, 396, 782–790.
58. N. D. Lepley and N. A. Holzwarth, Phys. Rev. B, 2015, 92, 214201.
59. S. Wenzel, S. J. Sedlmaier, C. Dietrich, W. G. Zeier and J. Janek, Solid State Ionics, 2018, 318, 102–112.
60. 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.
61. D. M. Anstine and O. Isayev, J. Phys. Chem. A, 2023, 127, 2417–2431.
62. J. Thomas, W. Baldwin, G. Csanyi and C. Ortner, Nonlinearity, 2025, 38, 095024.
63. J. L. Weber, R. D. Guha, G. Agarwal, Y. Wei, A. A. Fike, X. Xie, J. Stevenson, B. Santra, R. A. Friesner, K. Leswing, M. D. Halls, R. Abel and L. D. Jacobson, Efficient Long-Range Machine Learning Force Fields for Liquid and Materials Properties, arXiv, 2025, preprint, arxiv:2505.06462,  DOI:10.48550/arxiv.2505.06462.

---