Thermodynamics and Kinetics of the Cathode-Electrolyte Interface in All-Solid-State Li-S Batteries.

Lithium-sulfur batteries (LSBs) are among the most promising energy storage technologies due to the low cost and high abundance of S. However, the issue of polysulfide shuttling with its corresponding capacity fading is a major impediment to its commercialization. Replacing traditional liquid electrolytes with solid-state electrolytes (SEs) is a potential solution. Here, we present a comprehensive study of the thermodynamics and kinetics of the cathode-electrolyte interface in all-solid-state LSBs using density functional theory based calculations and a machine learning interatomic potential. We find that among the major solid electrolyte chemistries (oxides, sulfides, nitrides, and halides), sulfide SEs are generally predicted to be the most stable against the S8 cathode, while the other SE chemistries are predicted to be highly electrochemically unstable. If the use of other SE chemistries is desired for other reasons, several binary and ternary sulfides (e.g., LiAlS2, Sc2S3, Y2S3) are predicted to be excellent buffer layers. Finally, an accurate moment tensor potential to study the S8|β-Li3PS4 interface was developed using an active learning approach. Molecular dynamics (MD) simulations of large interface models (>1000s atoms) revealed that the most stable Li3PS4(100) surface tends to form interfaces with S8 with 2D channels and lower activation barriers for Li diffusion. These results provide critical new insights into the cathode-electrolyte interface design for next-generation all-solid-state LSBs.

Introduction

Lithium–sulfur batteries (LSBs) have garnered immense interest as a potential alternative to lithium-ion batteries (LIBs) for future energy storage needs due to their extremely high theoretical capacity (1672 mAh g–1)[1,2] and low cost due to the use of abundant sulfur.[3] However, commercialization of LSBs still requires solving several challenges that are intrinsic to this technology. While other challenges such as the electrically insulating nature and large volume change during cycling (∼80%) of the S cathode are being addressed via structure and composition engineering,[4−6] the polysulfide shuttling problem, which is considered to be the most detrimental factor to the cell performance, has not been completely mitigated.

Polysulfide shuttling refers to the concentration gradient-driven shuttling of soluble higher order polysulfides and insoluble lower order polysulfides between the cathode and anode through the liquid electrolyte medium, resulting in irreversible loss of the cathode and anode material. While the formation of polysulfides can be mitigated to some extent with alternative cathode chemistries, these solutions typically result in loss of energy density.[7−9] Replacing the liquid electrolyte with a solid electrolyte (SE) has been proposed as a potential solution to address the polysulfide shutting problem in LSBs.[10] For example, gel polymer electrolytes use polar molecules such as polyvinylidene fluoride[11] (PVDF) and poly(methyl methacrylate)[12] (PMMA) to trap dissolved polysulfides. In solid polymer electrolytes, the polysulfide migration is inhibited by their adsorption onto inorganic filler surfaces such as Al2O3[13] and TiO2.[14] Finally, composite SEs made by dispersing inorganic SEs in poly(ethylene oxide) (PEO)-based electrolytes[15−17] have been found to have a better ability to inhibit polysulfide movement and improve interfacial stability and mechanical properties when compared to single all-solid polymer or inorganic SEs.[10] In addition, an all-solid-state LSB architecture may yield other potential advantages of improved system-level energy density and safety.[18,19]

In contrast to all-solid-state LIBs,[20−27] the interfaces in all-solid-state LSBs, especially that between the SE and the S cathode, have not been extensively studied. The volatility of sulfur makes experimental characterization of these interfaces using traditional tools such as scanning electron microscopy and transmission electron microscopy challenging.[28−30] Nevertheless, existing studies suggest that these interfaces might not be stable. Zheng et al.[31] found that the Li6PS5Cl/S interface undergoes reductive degradation to form Li3P, LiP, Li3P7, or LiP7 along with Li2S and LiCl. It also undergoes oxidative degradation, resulting in thiosulfate oligomerization, which releases polysulfides that lead to sluggish Li-ion transport and cathode overpotential. An in situ study of a PEO-Li6.75La3Zr1.75Ta0.25O12 polymer ceramic composite interface with sulfur also showed the formation of several polysulfides, which migrate across the SE toward the Li anode.[32] While ab initio molecular dynamics (AIMD) simulations have been successfully used to probe the interfacial chemistries of solid-state (SS) interfaces in LIBs and LSBs,[33,34] their high computational cost severely restricts the size and the time scale of interfacial simulations. For instance, Camacho-Forero and Balbuena[34] found that the β-Li3PS4/cathode interface exhibits oligomerization of PS43– along with formation of various polysulfides within a 50 ps time scale, but interfacial strains of >10% were used to reduce the model size.

In this work, we aim to comprehensively investigate the electrochemical and chemical stability of cathode/SE interfaces in LSBs. We have intentionally not included the Li anode in our study given that this interface, which is no different from that in LIBs, has been extensively studied in prior works.[20−23,33] Using thermodynamic approximations,[35] we evaluated a large number of SEs with diverse anion chemistries (from oxides to sulfides to halides) for their interfacial stability against the charged S and discharged Li2S cathode. We find that sulfide SEs have the best (electro)chemical stability against the cathode in LSBs among the different classes of SEs. If the use of oxide/halide SEs in LSBs is desired for other reasons (e.g., potential stability against a Li metal anode), we also identified several sulfides that can potentially be used as buffer layers between the cathode and oxide/halide SEs. Finally, we developed a highly accurate machine learning interatomic potential (ML-IAP) via active learning for the β-Li3PS4|α-S8 model interfacial system. Molecular dynamics (MD) simulations on realistic large-scale (thousands of atoms) interfacial simulations using the developed moment tensor potential (MTP) reveal the formation of LiS and S species at the interface. The reaction products were in general found to increase the activation energy (Ea) of Li-ion diffusion at the interface, but the type of β-Li3PS4 surface interfaced was found to be the predominant factor that determines the interfacial Li-migration barrier.

Methods

Material Selection

Table presents the list of all the SEs and buffer layer materials studied in this work. Where available, the crystal structures and precomputed energies were obtained from the Materials Project (MP).[36] Otherwise, the initial crystal structures were obtained from the Inorganic Crystal Structure Database (ICSD).[37]Table S1 provides the sources of SEs considered in this study. We have attempted to be comprehensive in the selection of SE materials, covering all major anion chemistries. For the cathode, we considered both the fully charged (α-S8, mp-96) and discharged state (Li2S, mp-1153). For buffer layers, we considered all major experimentally reported oxide buffer layers[38] along with their sulfide analogues with an MP energy above convex hull (Ehull) of <30 meV. Tables S2 and S3 provide the sources and Ehull values of the oxide and sulfide buffer layers considered in this study.

Table 1

List of Solid Electrolytes and Buffer Layers Studied in This Worka

category	anion chemistry	compounds
solid electrolytes	oxide	Li1.3Al0.3Ti1.7(PO4)3,[39] Li7La3Zr2O12 (LLZO), Li0.33La0.56TiO3 (LLTO),[40] Li3.5Zn0.25GeO4,[41] Li2PO2N[42]
	nitride	Li3N
	sulfide	Li4GeS4, Li10GeP2S12 (LGPS), Li3PS4 (LPS),[43] Li7P3S11, Li6PS5Cl (LPSCl)[44]
	halide	Li3YCl6 (LYC), Li3YBr6 (LYB)
buffer layers	oxide	Li4Ti5O12, LiNbO3 (LNO), LiTaO3 (LTO), Li2Ti2O5, Li4Ti5O12, Li2ZrO3, Li2SiO3, Li3PO4, Li4TiO4, Li2TiO3, Li8Nb2O9, Li3NbO4, LiNb3O8, Li8SiO6, Li4SiO4, Li2Si2O5, Li5TaO5, Li3TaO4, LiTa3O8, Li4P2O7, LiAlO2, Li3BO3, LiH2PO4, LiTi2(PO4)3,LiBa(B3O5)3, LiPO3, LiLa(PO3)4, LiCs(PO3)2, Al2O3, ZnO, CdO, Sc2O3, Y2O3, La2O3, SiO2, TiO2, ZrO2, HfO2, Nb2O5, Ta2O5
	sulfide	Li2SiS3, Li3PS4, Li4TiS4, Li3NbS4, LiAlS2, Li3BS3, Al2S3, ZnS, CdS, Sc2S3, Y2S3, La2S3, SiS2, TiS2, ZrS2, HfS2

The corresponding material identifiers are provided in Tables S1, S2 and S3.

For phase diagram construction, the precomputed energies from MP[36] were used where available. In cases where the material was not part of the MP database (mostly SEs), the structures were obtained from the ICSD database, and DFT computations were carried out. For disordered structures, structure enumeration was performed on supercells to obtain ∼100 ordered candidates, and the lowest energy structure was used for further analysis.

Thermodynamic Approximations to Interfacial Stability

The electrochemical stability and chemical stability between SEs and LSB cathodes in charged (α-S8) and discharged state (Li2S) were studied using the well-established fast diffusion of Li-ion and multispecies equilibrium approximations, respectively.[35,45−50] Here, a brief summary is provided, and interested readers are referred to those previous works for details.

To estimate chemical stability, the reaction energy is calculated by constructing a pseudobinary phase diagram between the SE and cathode and determining the ratio that results in the most negative reaction energy. The assumption here is that the two materials react to form the most favorable products under thermodynamic equilibrium. The reaction energy is given byHere cSE and ccathode are phases of the SE and cathode, which are the terminal compounds of the pseudobinary phase diagram. E(cSE) and E(ccathode) are their respective density functional theory (DFT) energies. Eeq(c) is energy of phase equilibria at a given composition c, and N is the number of atoms involved in the reaction used as the normalization factor. ΔE(cSE, ccathode) is the estimate for thermodynamic (equilibrium) stability of the SE and cathode.

The electrochemical stability of the SE/cathode interface was evaluated using the grand potential phase diagram with respect to Li. The interface can be approximated as the SE composition in contact with a Li sink or source. For a system open with respect to Li, the relevant thermodynamic potential is given by the grand potential Φeq, as follows:

The change in the grand potential (ΔΦ) of the SE/cathode interface is given by eq and provides an the estimate of the thermodynamic stability under cycling conditions.

It should be noted that the normalization factor Ngc in eq , unlike N in eq , is the total number of atoms excluding Li. The estimated percentage volume changes (ΔV) for thermodynamic reactions were computed using the DFT relaxed volumes of the reactants and products and weighting them in accordance with the ratios in the balanced reaction equation. A positive ΔV indicates a net increase in the volume due to the reaction and suggests the possibility of stress buildup and cracking of active material at the interface. A negative ΔV indicates a net decrease in the volume due to the reaction, possibly leading to the formation of voids at the interface.

DFT Calculations

All DFT calculations were performed using the Vienna Ab initio Simulation Package (VASP) with the projector augmented wave (PAW) approach.[51,52] For thermodynamic analysis, the bulk structure relaxations were done using the Perdew–Burke–Ernzerhof (PBE) generalized gradient approximation (GGA) functional,[53] and the calculation input parameters were kept consistent with MP. The calculations were spin polarized with an energy cutoff of 520 eV and k-point density of 64/Å–3.

For interface construction, the bulk structures were relaxed with stricter energy and force convergence criteria of 10–5 eV and 0.02 eV/Å, respectively. Additionally, for systems containing α-S8 (bulk and interface), the relaxation was performed using the optB88 van der Waals (vdW) functional.[54] The energy cutoff was reduced to 400 eV, and a looser energy convergence criterion of 10–4 eV was used.

All the DFT, AIMD simulations, and output analysis were performed using automated workflows[55] developed using Python Materials Genomics (Pymatgen)[56] and Fireworks[57] packages.

Interfacial Moment Tensor Potential Development

In recent years, machine learning (ML) of the potential energy surface as a function of local environment descriptors[58−65] has emerged as an automatable approach to develop interatomic potentials (IAPs) for complex chemistries, including SE materials,[66,67] with near DFT accuracy. In this work, we have adopted the moment tensor potential (MTP) formalism,[61] which has been shown to provide a good balance between accuracy and computational cost.[68] Interested readers are referred to those works for details on the MTP formalism.[61]

The α-S8|β-Li3PS4 interface was selected as our model interface of interest. To capture the full complexity of local environments and bonding across the different phases of interest (Li2S, α-S, β-Li3PS4) and at the interface, we adopted a sequential passive learning and active learning workflow adapted from Novikov et al.[69] to fit the MTP, as shown in Figure . In the passive learning stage, training structures were generated by performing non-spin-polarized AIMD simulations using the NVT ensemble on supercells of relaxed bulk structures of β-Li3PS4, cubic Li2S, and α-S8 with the PBE functional. Supercells with cell parameters of at least 10 Å were used to minimize the effects of periodic images. A Γ-centered 1 × 1 × 1 k-point mesh was used with a plane wave energy cutoff of 400 eV. Temperature control was achieved using a Nose–Hoover thermostat,[70,71] and the time step was fixed at 2 fs. During the simulation, the 0 K relaxed structures were heated to the target temperature with a temperature gradient of 0.5 K/fs followed by an equilibration period of 10 ps and a production run of 20 ps. The simulations were carried out at 300, 600, 900, 1200, and 2000 K, and 150 snapshots were collected from each temperature at an interval of 0.2 ps, totaling 2250 structures. Further, 10 snapshots were collected at an interval of ∼3 ps at each temperature for Li2S and β-Li3PS4 and volumetrically strained (−10 to +10%, 2% interval) to generate an additional 1000 structures. The training set was further augmented by adding ∼100 slab structures per compound, which were generated by a similar straining procedure of stoichiometric symmetric slabs of these compounds.

Figure 1

Passive and active learning flowchart for interfacial MTP development.

Static self-consistent calculations were then performed using the optB88 vdW functional[54] on the training structures to obtain accurate energies, forces, and stresses. The optB88 vdW functional was selected for the energy and force evaluations given that we expect that vdW forces to play a substantial role in these materials.[66] Further, the lattice parameters calculated with the optB88 functional for Li2S, S, and Li3PS4 are in much closer agreement with experimental values (see Table S4). The two key parameters that control the trade-off between accuracy and computational cost of MTPs are the radius cutoff (Rcutoff) and the maximum level (levmax). MTPs with different radius cutoffs (Rcutoff) and maximum levels (levmax) were trained and selected based on convergence of mean absolute error (MAE) for energies (MAEenergies) and forces (MAEforces) of test structures with respect to levmax. The training set consisted of 384 training structures, which are tabulated in Table S5. For MTP fitting, a 90:10 training:test split was used with the weights of 100:1:0.1 for energy, force, and stresses. Three MTPs with an Rcutoff of 5, 6, and 7 Å with a levmax of 20 were chosen for improvement with active learning (AL) (see Figure S1).

For the AL stage, interface models were constructed with slab structures of the different surfaces of β-Li3PS4 and α-S8 using the algorithm proposed by Lazić.[72]Figure S2 shows the slab unit cells used for interface cell construction. The values of DFT surface energies are tabulated in Table S6. For DFT calculations and MTP model training and validation, smaller interfacial models (<250 atoms) were constructed by allowing interfacial mean absolute strains of up to 10% and cell parameters of at least 10 Å along the direction of the interface. Figure S3 shows the seven distinct α-S8|β-Li3PS4 interfaces that satisfy these criteria. Table shows the surfaces used to construct the interfaces, their DFT-computed interfacial energies, and the reference names. The complete details of the parameters associated with the interfaces and their construction are tabulated in Table S7.

Table 2

Interfaces Constructed Using Different Surfaces of β-Li3PS4 and α-S8 along with DFT Interfacial Energies

β-Li3PS4 surface	α-S8 surface	interface energy (J/m2)	name
100	001	33.16	S8(001)|Li3PS4(100)
	111	22.88	S8(111)|Li3PS4(100)
001	001	34.15	S8(001)|Li3PS4(001)
	111	41.89	S8(111)|Li3PS4(001)
010	001	39.12	S8(001)|Li3PS4(010)
	001	40.82	S8(001)|Li3PS4(010)(2)
	111	57.58	S8(111)|Li3PS4(010)

A detailed description of the AL process is provided in ref (69). Here, MD simulations were performed at four temperatures (300, 500, 700, and 900 K) on the seven interfacial structures, with three levels of isotropic strain (−5%, 0%, +5%) applied, resulting in 4 × 7 × 3 = 84 MD simulations. The extrapolation grade (γ) is used to select the distinct “new” structures that augment the training set. γ is a measure of the degree of extrapolation the MTP performs to estimate the energy of a structure[73] and serves as a proxy for energy and force estimation errors.[69] In practice, two threshold values of γ, γselect and γbreak, are defined. Structures with γselect < γ < γbreak are considered reliable extrapolation and are added to a preselection pool, whereas structures with γ > γbreak are deemed to be “risky” and result in termination of simulation.[69] The γselect and γbreak were initially set to 2 and 10, respectively. The γbreak value was dynamically changed during iterations between 10 and 2 to ensure no more than 40 000 configurations were selected in any AL iteration. The MD simulations that had no contribution to the preselected structure pool of an iteration were discontinued in future iterations. In addition to reducing computational cost, this tailors the AL workflow toward capturing more diverse and complex local environments, i.e., those generated at higher temperatures and reactive interfaces. Static calculations were performed on the selected structures (a subset of the preselected structure pool) to update the training set and obtain the next iteration of MTP. This process was repeated until the resulting MTP met four conditions: (i) 100% MD reliability (no termination within 1 ns); (ii) fewer than 50 structures were added to the selection pool; (iii) test MAEenergies < 5 meV/atom; and (iv) test MAEforces < 0.2 eV/Å. The training and evaluation of MTPs were performed using the MLIP[61] and Materials Machine Learning[74] (maml) packages, and MD simulations were performed using the LAMMPS[75] package.

Molecular Dynamics Simulations and Analysis

To study the kinetic stability of SE/cathode interfaces, MD NPT simulations were performed using an accurate MTP on large interface models (>1000 atoms) at a 5 ns time scale. An elevated temperature of 600 K was used to speed up the possible chemical reactions at the interface. A 2 fs time step was used, and the interfaces were heated and equilibrated for 40 ps each before the production run of 5 ns. The interface structure evolution and reactions were analyzed using the clustering algorithm implemented in OVITO visualization software.[76] The molecular nature of α-S8 and the interface reaction products allows for the use of a bond connectivity based clustering algorithm to decompose the interfacial structure into molecular clusters.

A cluster is defined as a set of connected points within which every atom is within the reach of every other atom through at least one continuous path. The path is made of atoms and their bonds with the adjacent atoms. Element pairwise bond length cutoff criteria were used to consider a pair of atoms bonded. The cutoffs were set at values 10% higher than the equilibrium bond lengths to account for bond expansion due to heating of the structure in MD simulations. Figure S4 provides examples of clustering of atoms based on the bond length cutoff criteria.

For Li diffusion analysis, the diffusivity (D*) was calculated by performing linear fitting of mean square displacement (MSD) of a Li ion with time using the Einstein relation[77]Here, N is the total number of Li ions and Δr(t) is displacement of the ith ion at time t. For calculation of surface/interfacial diffusivities, only Li atoms with both start and end positions of hopping inside the surface/interfacial region were considered for MSD calculation.[78−81] The Arrhenius plot was then obtained to calculate the bulk/surface/interface activation barrier, Ea.

Results and Discussion

Thermodynamic Analysis of SE/Buffer/Cathode Interfaces in LSBs

Chemical Stability of SE/Cathode Interfaces

Figure shows the reaction energy (left) and corresponding volume change (right) for the different SE/electrode under consideration. The complete list of reaction products of all SE/cathode pairs is provided in Table S8. Li3N/S8 is found to be the most unstable SE/electrode pair with a very high reaction energy of −900 meV/atom. This results in the formation of Li2S and the evolution of N2, which might present a safety hazard. All oxide SEs are found to be unstable against S8 with similar reaction energies and a volume shrinkage of ∼250–300 meV/atom and ∼8–10%, respectively. The reaction products generally include metal sulfides such as LaS2, ZrS3 (TiS3), and Li2S. Among the halide SEs, the LYC/S8 interface is predicted to be chemically stable. LYB reacts with a small reaction energy of 22 meV/atom and a negligible volume increase of ∼1% to form binary compounds such as Y2S3, LiBr, and SBr.

Figure 2

Reaction energy (left) and percentage volume change (right) due to reaction for different electrolyte/electrode pairs.

The sulfide SE/S8 interfaces show the best chemical compatibility except in the case of Li7P3S11, which reacts to form Li3PS4 and P2S7. The discharged Li2S cathode is generally more chemically stable against all oxide and sulfide SEs compared to the charged S8 cathode. Surprisingly, the discharged Li2S cathode is chemically less stable against halide SEs when compared to the charged S8 cathode. This is due to the formation of highly stable compounds such as LiX (X = Cl, Br) and LiYS2.

Electrochemical Stability of SE/Cathode Interfaces

Figure shows the electrochemical reaction energy (left) and corresponding volume change (right) for SE/α-S8 interfaces between 0 and 2.25 V. In general, the qualitative observations about the relative chemical stabilities of different SE anion chemistries also apply to the relative electrochemical stabilities. The Li3N/S8 interface is electrochemically stable up to 2 V, after which the interface is predicted to decompose to Li and S7N with a high reaction energy of −865 meV/atom and volume change of about 2%. The halide SEs interfaces show very limited electrochemical stability and start to react at 0.25 V to form LiX (where X = Cl, Br) and YS with a high reaction energy of >500 meV/atom. However, the volume changes are <10%. Similarly, the oxide LLZO and LLTO interfaces are electrochemically unstable at all voltages. Li1.3Al0.3Ti1.7(PO4)3 and Li3.5Zn0.25GeO4 have slightly larger stability windows of 0–1.25 V and 0–1 V, respectively. Finally, the sulfide SE interfaces are electrochemically stable at all voltages with the exception of Li7P3S11 at 1.5 V, where it reacts with S8 to form Li3PS4.

Figure 3

Electrochemical reaction energy (left) and corresponding volume change (right) for α-S8.

Compared to the charged S8 cathode, the discharged Li2S cathode/SE interfaces have identical reactions with SEs, but the reaction energies and resultant volume changes are in general predicted to be significantly lower than α-S8, as shown in Figure S5.

Chemical Stability of SE/Buffer Layer Interfaces

In general, oxide and halide SEs are predicted to be both chemically and electrochemically unstable against LSB cathodes to varying degrees. Nevertheless, there may be potential situations where the use of an oxide or halide SE may be desired. For instance, oxide SEs may potentially provide better chemical and mechanical compatibility with the Li metal anode. Here, we investigate potential oxide or sulfide buffer layers that can be used to stabilize interfaces. Figures and 5 show the predicted reaction energies and corresponding percentage volume changes for SEs and cathode interfaces with oxide and sulfide buffer materials.

Figure 4

Reaction energy (top) and volume change (bottom) for different electrolyte and electrode with oxide buffer layer material pairs.

Figure 5

Reaction energy (left) and volume change (right) for different electrolytes and electrodes with sulfide buffer layer material.

As can been seen, like oxide SEs, all the oxide buffer materials are predicted to be incompatible with a S8 cathode. However, some sulfides show promise of being used as buffer layers. For example, some binary and ternary sulfides such as Sc2S3, Y2S3, Li3PS4, and LiAlS2 are predicted to have chemical compatibility with both Li2PO2N and the charged and discharged LSB cathodes. The binary sulfides also exhibit good chemical compatibility with halide SEs. Li3PS4 can also be used to buffer the Li1.3Al0.3Ti1.7(PO4)3 interfaces as well.

Atomistic Studies of the S8/β-Li3PS4 Interface

Given that sulfide SEs appear to be the most stable against the S8 cathode in LSBs, we have selected the β-Li3PS4/S8 interface as our model system for in-depth studies of interfacial kinetics via MD simulations using a fitted MTP. In the following sections, we will first provide a validation of our fitted MTP before analyzing the results of the MD simulations.

MTP Validation

Figure a and b show the evolution of the test MAEenergies and MAEforces for different interfaces during the AL process. The test structures were generated by production MD simulations at 300 and 600 K using the MTP trained at the last (25th) iteration of AL. Test structures were collected at 2 ps intervals during the heating and equilibration stages and at 20 ps intervals during the production stage. It may be observed that the test MAEenergies and MAEforces converge to below 5 meV atom–1 and 150 meV Å–1, respectively, for all interfaces after about 15 AL iterations. The fitted MTP is also able to reproduce the lattice parameters and relative energy differences for the α, β, and γ polymorphs of Li3PS4, as shown in Table S9.

Figure 6

(a) MAEenergies and (b) MAEforces for test structures obtained from different interfaces during active learning iterations. The test structures were collected from MD simulations performed at 300, 500, 700, and 900 K of the last iteration AL to ensure no structures were part of the training set (added back during AL).

Figure a–c compare the test MAEenergies and test MAEforces as well as calculated MAEinterface energies for all the interfaces. Across all interfaces, the MAEenergies and MAEforces are consistent with the results obtained during the AL iterations. The MAEinterface energies are within 0.15 J m–2. It is important to note that the MTP was not trained on MD structures from 600 K during AL, and yet the errors for 600 K structures are low. This indicates that the MTP is unlikely to be overfitted.

Figure 7

MAE of MTP-predicted (a) energies, (b) forces, and (c) interfacial energies relative to DFT values for test structures collected from 300 and 600 K production run MD simulations.

Table compares the DFT and MTP-computed cell constants, mechanical properties, and surface energies for β-Li3PS4 and α-S8. In general, we observe that the errors in the cell constants and mechanical properties are extremely low. However, the MTP tends to underestimate surface energies for both β-Li3PS4 and α-S8. The MTP-predicted surface energies for α-S8 are well within 0.05 J/m2, possibly due to their lower bonding complexity as well as the lower absolute values of the surface energies. The errors in the MTP-predicted surface energies of the more complex and higher-energy β-Li3PS4 surfaces are somewhat higher, but in line with the MAE in interface energies during the AL iterations.

Table 3

Comparison of DFT- and MTP-Computed Cell Parameters, Mechanical Properties, and Surface Energies of β-Li3PS4 and α-S8a

	β-Li3PS4		α-S8
quantity	DFT	MTP	DFT	MTP
Cell Parameters
a (Å)	13.07	12.99 (−0.08)	10.33	10.31 (−0.02)
b (Å)	8.13	8.09 (−0.05)	12.83	12.80 (−0.03)
c (Å)	6.26	6.23 (−0.03)	24.5	24.44 (−0.06)
Mechanical Properties
bulk modulus (GPa)	22	29.61 (7.61)	14.34	20.25 (5.91)
shear modulus (GPa)	11	14.81 (3.81)	7.28	11.77 (4.48)
Poisson’s ratio (GPa)	0.27	0.32 (0.05)	0.28	0.28 (0.00)
Surface Energies (J/m2)
surface index
100	0.371	0.321 (−0.050)	0.172	0.151 (−0.021)
001	0.608	0.449 (−0.159)	0.188	0.157 (−0.031)
010	0.620	0.560 (−0.060)	0.185	0.172 (−0.013)
111			0.153	0.132 (−0.021)

Values in parentheses show the absolute error of MTP-computed values relative to the DFT values.

MD Simulations of β-Li3PS4/S8 Interfaces

Figure a–f show the time evolution of cluster size (nc) distribution of β-Li3PS4 interfaces with the S8 surfaces with the highest (left panels) and lowest (right panels) surface energies, i.e., the (001) and (111) surfaces, respectively. It should be noted that no phase transformation from β-Li3PS4 to either γ-Li3PS4 or α-Li3PS4 was observed within the MD simulation time scale, as such transformations require a kinetically activated rearrangement of PS4 and LiS4/LiS6 polyhedra. In general, it can be seen that the percentage of clusters with nc < 8 increases steadily over the simulation time with a concurrent decrease in the percentage of clusters with nc = 8 (S8 molecules), indicating the breakdown of S8 rings into smaller S rings/chains. However, there is a noticeable difference in the evolution of α-S8 (001) and (111) based interfaces. The clusters with nc ≠ 8 form sooner and at much higher concentrations at the S(001)-based interfaces (Figure a–c). There is also a significantly smaller percentage of clusters with nc > 9 at the S(111)-based interfaces (Figure d–f). Figure a–f show the average cluster size n̅c distribution and composition within each cluster size for the β-Li3PS4 interfaces with the S8 surfaces with the highest and lowest surface energies, i.e., the (001) and (111) surfaces, respectively. The averaging was done over the last 1 ns of simulation. The molecular structures constituting a given cluster size are shown on top of each cluster size in the same order as their percentage subcomposition shown in cluster size bars. For example, in Figure a, the nc = 3 clusters consist of S3 and LiS2 type molecules at ∼1% and ∼0.1% composition, respectively. The clusters that were present at less than 0.1% concentration are classified as “others”. It is also important to note that Li atoms remain bound to β-Li3PS4 and any molecular cluster containing Li should be thought of as S rings/chain/atom interacting with the surface Li atoms of β-Li3PS4 and not as free-floating Li-polysulfides.

Figure 8

Time evolution of cluster size distribution in (a) S8(001)|Li3PS4(001), (b) S8(001)|Li3PS4(010), (c) S8(001)|Li3PS4(100), (d) S8(111)|Li3PS4(001), (e) S8(111)|Li3PS4(010), and (f) S8(111)|Li3PS4(100) interfaces.

Figure 9

Averaged cluster size distribution and composition at (a) S8(001)|Li3PS4(001), (b) S8(001)|Li3PS4(010), (c) S8(001)|Li3PS4(100), (d) S8(111)|Li3PS4(001), (e) S8(111)|Li3PS4(010), and (f) S8(111)|Li3PS4(100) interfaces.

While there are some differences in the percentage composition of the different interfaces, the qualitative trends in cluster distribution profiles are remarkably similar. Each cluster size predominantly consists of S rings/chains at 2–3% followed by the presence of a small percentage (<1%) of Li polysulfide-like (Li–S) molecules. The Li–S molecules are formed by the interaction of S rings/chains with the Li ions from β-Li3PS4 present at the interfaces.

Based on the bond lengths and bond angles (accounting for thermal expansion of bonds), the nc = 3 clusters strongly resemble the known structures of S3 (S3 trimer) and solid-state c-Li2S. The nc = 2 S2 clusters also closely resemble the known S2 dimer structure. While the presence of shorter Li–S (x < 8) molecules can be expected at the interface, the presence of larger Li–S (x > 8) molecules at a similar percentage composition is interesting. These larger Li–S (x > 8) molecules are formed by the fusion of smaller S chains followed by absorption onto the β-Li3PS4 surface at the interface. This might indicate a fundamental difference in the interface composition between liquid and solid-state Li–S batteries.

From Figure d–f, it may be observed that the percentage of smaller clusters (nc < 8) is much lower (<1%) for the S8(111) interfaces compared to the S8(001) interfaces; that is, the S8(111) interfaces are less reactive than the S8(001) interfaces. This is likely due to the fact that the (111) surface is significantly lower in energy and, hence, more stable than the (001) surface. This is also supported by the plateauing of the percentage of cluster sizes at ∼2% in Figure S7a–f, which is not observed in the case of the (001) surface (Figure S6a–f). As found in the time evolution analysis, there is a marked decrease in the presence of larger Li–S (x > 10) polysulfide resembling structures. However, the types of shorter S rings/chains and Li–S (x < 8) polysulfide structures are similar to that of α-S8(001) interfaces. It is interesting to note the similarities and differences in the interfacial reaction products from a previous AIMD study on the same kind of interface.[34] While the formation of a few smaller S species such as S dimers and trimers was shown, this study found no evidence of the formation of larger S and LiS types of species. This difference, among other factors, such as interface morphology and interfacial strain, might be due to the limited system size and time scale of the simulations possible using ab initio methods. The larger system size and time scale of MTP-driven MD simulations allow for extended periods of chemical interactions of various species at the interface, thus allowing for the formation of larger polysulfide species. This also eliminates the need for using less realistic interfaces with high interfacial strain (used to limit system size), which also presumably impacts the interfacial reactivity.

Figure a and c show the time evolution of the spatial distribution of non-S8 and S8 clusters, respectively, as a function of distance along the direction perpendicular to the interface for the S8(001)|Li3PS4(001) interface. The interface structure is shown in Figure b. It can be seen that the thickness of the reaction zone (interphase) increases initially, after which it plateaus. There is a concurrent decrease in the concentration of S8 ring atoms in the reaction zone. The bulk-like regions of α-S8 are unaffected. It is important to note that while the interphase region grows, there is no evidence of diffusion of reaction products (Li polysulfides) into the Li3PS4 bulk region within the simulated time scale. This is true for all the interfaces under consideration, as evident from Figures S8a–f. Thus, crystalline β-Li3PS4 is effective in blocking polysulfide diffusion. It should be noted that, unlike a full cell battery, there is no applied electrochemical driving force for the Li from β-Li3PS4 to diffuse into α-S8 except for the difference in the Li chemical potential between the compounds. Further, the presence of grain boundaries at the interface will also impact the interfacial reactivity, interphase thickness, and polysulfide mobility, which will be the focus of a future study.

Figure 10

(a) Time evolution of normalized atom count of interface reaction products (cluster size ≠ 8) as a function of the distance for the S8(001)|Li3PS4(001) interface. The center of the bulk Li3PS4 region is the zero reference. (b) S8(001)|Li3PS4(001) interface cell with highlighted interface region. The Li, P, and S atoms are represented by green, purple, and yellow atoms, respectively. (c) Time evolution of normalized atom count of S8 atoms as a function of the distance from the S8(001)|Li3PS4(001) interface. (d) Time evolution of the interphase layer thickness in α-S8|β-Li3PS4 interfaces. The average interphase thickness is averaged over the last 2 ns.

Figure d shows the evolution of interphase thickness for the interfaces during the simulation time scale. The interphase thickness is averaged over the two interfaces present in the interface cell. The time-averaged interphase thickness is calculated over the last 2 ns of the simulation. It can be seen that the interphase grows very rapidly within the first 100 ps for all the interfaces, after which there are only minor fluctuations. The α-S8(001) surface based interfaces with higher reactivity form a considerably thicker interphase layer of 6.5–8 Å. The more stable α-S8(111)-based interfaces have an interphase thickness of 3.5–4.5 Å. Hence, surface stability can have a profound effect on the interfacial reactivity and, correspondingly, the thickness of the interphase layer. These interphase thicknesses can be considered as the spatial limits of interface-related phenomenon and are used as the “interface region” for the analysis of ionic diffusivities in this region.

Figure a,b show the Arrhenius plots and evolution of the interfacial areal number density of Li nLi′ during simulation time scale for the interface region of the different α-S8|β-Li3PS4 interfaces. For comparison, we have computed the Arrhenius plots of bulk and different surfaces of β-Li3PS4 in Figure S9a. The Ea of bulk β-Li3PS4 is estimated to be 0.338 eV, in agreement with the previously reported value.[82−84] The (001), (010), and (100) surfaces of β-Li3PS4 are estimated to have the Ea values of 0.373, 0.290, and 0.320 eV, respectively. The surface Ea was obtained from the diffusivity of surface Li atoms of β-Li3PS4 slabs present in the same interphase thickness range as their derived interfaces. Figure S9b shows the evolution of nLi′ for different β-Li3PS4 surfaces during the simulation time scale. The β-Li3PS4 surfaces with lower surface Ea compared to bulk, i.e., (010) and (100) surfaces, have an estimated nLi′ of ∼13.5 nm–2, which is almost twice that of the higher Ea (001) surface with an nLi′ of 7.5 nm–2. The Li diffusion pathways for lower Ea surfaces are found to be two-dimensional, but those for the higher Ea (001) surface are predominantly one-dimensional (Figure S9c–e). Hence, it is inferred that surfaces with lower nLi′ and localized Li diffusion pathways tend to have higher surface Ea. In general, we find that the interphase formation results in an increase of Ea when compared to the parent β-Li3PS4 surface. However, this is also dependent on the relative Li enrichment/depletion of the interface. For example, S8(111)|Li3PS4(001) has the highest estimated interface Ea of 0.456 eV, which is ∼0.08 eV higher than the parent β-Li3PS4(001) surface. Additionally there is a marginal drop in the interface nLi′ value to 6.89 nm–2 from 7.5 nm–2. Li trajectories of this interface (Figure c) also reveal a disruption of the one-dimensional Li diffusion pathways. In contrast, the S8(001)|Li3PS4(001) interface containing the same β-Li3PS4(001) parent surface shows marginal increase in the Ea despite the formation of a thicker interphase. It can be observed that the nLi′ for this interface increases over time and averages ∼2 nm–2 higher than the parent surface. Similarly, the S8(001)|Li3PS4(100) interface exhibits a marginal drop in the Ea of ∼0.03 eV compared to the parent β-Li3PS4(100) surface and has the highest nLi′, which is a marginal increase of ∼1.5 nm–2. This is reflected in the marginal increase in the Li trajectories’ coverage when compared to the parent surface (Figure S9e), as seen in Figure d.

Figure 11

(a) Arrhenius plot and (b) time evolution of Li areal density at the interfacial regions in α-S8|β-Li3PS4 interfaces. The activation barrier for Li-ion migration and average Li areal density nLi′ at the interface region are provided. The Arrhenius plot for bulk β-Li3PS4 is also provided for reference. Li trajectories (colored in red) from MD simulations for (c) S8(111)|Li3PS4(001) and (d) S8(001)|Li3PS4(100) interfaces viewed along the c crystallographic axis (direction perpendicular to interface). Surface morphology of (e) Li3PS4(001) and (f) Li3PS4(100) surfaces.

Finally, it should be noted that while the changes in interface nLi′ and interface Li trajectories can be used to explain the changes in interface Ea, the type of β-Li3PS4 surface predominantly determines the Ea at the interface. Surfaces with higher Ea tend to form interfaces with high Ea. These differences can be traced to the morphology and Li diffusion topology of different β-Li3PS4 surfaces (Figures e,f). The (001) surface is corrugated, which results in the formation of one-dimensional Li diffusion channels with few opportunities for interchannel Li hopping. In contrast, the (100) surface (and (010) surface) is highly planar, resulting in a multitude of diffusion pathways that achieve high coverage of the surface, resulting in lower Ea.

Conclusion

To conclude, this work provides a comprehensive thermodynamic and kinetic analysis of the cathode/SE interface in LSBs. In general, the charged S8 cathode was predicted to be much more reactive than the fully discharged Li2S cathode. Among the different anion chemistries, sulfide SEs were found to be the most thermodynamically stable against the LSB cathode. Sulfides are also highly effective buffer layers if the use of other SE anion chemistries is desired for other reasons. We expect these findings to be generalizable to any new SEs or buffer layers discovered for either LIBs or LSBs. Through nanosecond MD simulations using a highly accurate MTP, it was demonstrated that the reactivity at the Li3PS4/S8 interface is related to the relative surface energies of the composite α-S8 surfaces; a higher energy surface tends to be more reactivity than a lower energy surface. Finally, while the formation of an interphase tends to increase activation energies for Li diffusion, the surface morphology and Li diffusion topology of the Li3PS4 surface were found to be the critical determinants of interfacial Li diffusion dimensionality and barriers. The most stable Li3PS4(100) surface tends to form interfaces with S8 with 2D channels and lower activation barriers for Li diffusion.