Theoretical study on stability and ion transport property with halide doping of Na₃SbS₄ electrolyte for all-solid-state batteries

Abstract

All-solid-state Na ion battery (ASS-NIB) is a new class of battery which is a potential alternative to the conventional all-solid-state Li ion battery. Herein, we focused on Na₃SbS₄, a reported candidate solid electrolyte for ASS-NIB application, and performed a comprehensive theoretical study, primarily based on density functional theory calculations, to evaluate (electro)chemical stability, defect chemistry, and Na ion transport property. The calculated results reveal that when in contact with a layered cathode compound (e.g., NaCrO₂), sulfur in Na₃SbS₄ tends to migrate across the interface, leading to interface atomic rearrangement, interface disordering and/or decomposition. The material is also predicted to decompose under reductive voltage conditions (0 V vs. Na⁺/Na), in agreement with experiment. Kinetic modeling for stresses and electron density distribution showed that the interface between Na₃SbS₄ and Na metal anode under typical electrodeposition surface roughness would lead to dendrite initiation and growth. Aside from controlling Na vacancy concentration, halide doping at the S site was predicted by DFT molecular dynamics (MD) calculations to directly affect the Na⁺ ion activation energy. This can be ascribed to the size modulation of the Na site-to-site pathway bottleneck. Cl and Br halide dopants are both determined to be promising for conductivity optimization, with low DFT-MD Na⁺ ion activation energy (∼0.1 eV at 4% Na vacancy). Thermodynamic analysis shows the possible factors that can influence the final conductivity of Na₃SbS₄: secondary phases and intrinsic defects. Overall, our findings offer valuable insights for the rational design of solid electrolytes.

---

Introduction

Due to cost and accessibility issues associated with Li and the ever increasing desire of the society to tap more into clean/renewable energy sources that have less carbon footprint, many research studies now focus on alternative technologies besides Li-ion rechargeable batteries, one being based on the highly abundant Na element.^1–4 Na is especially appealing cost-wise such as in stationary applications wherein gravimetric and form factors are not as crucial (e.g., in back-up power systems, smart grids, etc.). Material components in Na and Li ion batteries share many similar chemistries. As a result, the former has benefited in terms of accelerated development cycles by leveraging on accumulated knowledge and expertise from the latter.^5–22

Similar to the Li ion battery, the increasing demand for high safety in electric/electronic applications has also motivated many to rethink the current design of a typical sodium ion battery (NIB). One approach that is seriously considered as a solution is the replacement of the combustible conventional organic-/liquid-based electrolyte with a solid electrolyte (SE), shifting the current paradigm towards an all-solid-state Na ion battery (ASS-NIB) concept. Some of the candidate SE materials that have been explored are Na-β′′-alumina,²³ NASICON NaZr₂(PO₄)₃ and Na_3.4Sc_0.4Zr_1.6(SiO₄)₂(PO₄),^24,25 borohydride Na₂B₁₂H₁₂,²⁶ anti-perovskite Na₃OBr,²⁷ and sulfides such as Na₁₁Sn₂PS₁₂,²⁸ Na₃PS₄,²⁹ Na_3−xPS_4−xCl_x,³⁰ Na₃PSe₄ (ref. 31) and Na₃SbS₄.^32–34 Specifically, Na₃SbS₄ (NSS), a recently reported one, has shown great promise in terms of ease in processing,³³ ionic conductivity (i.e., order of 10⁻³ S cm⁻¹),^32,34 and material processing (i.e., less-intensive conditions are needed for high-elasticity sulfide-type SEs, as compared to stiff oxide-type SEs). The good interface mechanical behavior of NSS when in contact with a cathode (NaCrO₂) under cycling operation has also been reported.³³ NSS has also been found to be more chemically stable than a predecessor electrolyte Na₃PS₄, in particular against toxic H₂S gas formation from H₂O contact.^29,32

NSS was noted to undergo chemical decomposition when in direct contact with a Na metal anode (without cell biasing), forming decomposition phases such as Na₂S and Na₃Sb based on XPS analysis.³⁵ Such anode-side reactions can have profound impact on battery cell performance. Another issue that is yet to be studied is the electrodeposition stability with Na–SEs against initiation of dendrites on the anode surface.

In terms of ionic conductivity (σ), NSS was reported to vary by a factor of 1–2 between its two known phases, the low-temperature tetragonal (t-NSS) phase and the high-temperature cubic phase (c-NSS), with comparison being 1.77 × 10⁻³ S cm⁻¹vs. 2.80 × 10⁻³ S cm⁻¹, respectively.^32–34,36 The ion conduction mechanism in both phases was presumed to be governed by Na vacancies. A DFT-based calculation using a dilute vacancy model with the single hopping mechanism has predicted an enormously low local energy barrier (<0.1 eV) along the Na diffusion channel.³⁷ Since then, σ improvement has been attempted in several studies such as by applying pressure (1.00 × 10⁻³ S cm⁻¹ at 0.9 GPa),³⁸ isovalent substitution at the S²⁻ site by Se²⁻ (0.85 × 10⁻³ S cm⁻¹),³⁹ and aliovalent substitutions at the Sb⁵⁺ site by Sn⁴⁺ (0.50 × 10⁻³ S cm⁻¹),⁴⁰ W⁶⁺ (10⁻² S cm⁻¹),^41,42 and Mo⁶⁺ (10⁻³ S cm⁻¹).⁴³ Meanwhile, aliovalently doping the S site with halide anions is another interesting strategy which can modulate the vacancy concentration and diffusion pathway. In particular, the ion conduction channel size can be largely affected by variation in the anion sublattice which largely determines the packing in crystalline solids. The remaining spaces excluding the anion sublattice are Na sites, non-Na cation sites (e.g., Sb⁵⁺ in NSS), and interstitial sites.

Given the high expectation for NSS as a key SE material for next-generation ASSBs, we thus aimed in this work to further investigate it on three key issues that are crucial for practical application: (i) compatibility with the cathode (i.e., crossing tendency of chemical species across the cathode–electrolyte interface that may lead to increased impedance), (ii) stability with the Na metal anode from the viewpoint of decomposition and Na dendrite growth, and (iii) ion transport with respect to halide anion substitution at the S site. Our results not only provide theoretically determined properties that can be compared directly with experimental values, but they can also be used when formulating possible material optimization strategies for solid electrolytes in general.

Computational details

Geometry optimization

DFT calculations were performed using the VASP code^44,45 which implements the projector augmented wave (PAW) formalism for the ion–electron interactions.^46–48 The generalized gradient approximation (GGA) as parameterized by Perdew, Burke, and Ernzerhof (PBE) was employed as the functional for the exchange–correlation energy.^49,50 Except for the Na pseudopotential which includes semi-core p states as valence states, standard pseudopotentials were used for Sb, S, and dopant anions F, Cl, Br, and I. The cutoff for kinetic energy was set to 520 eV and a k-point mesh of at least 1000 was assigned with the Monkhorst–Pack grid scheme.⁵¹ All calculations were carried out under spin-polarized conditions. Initial structural data were taken from the Inorganic Crystal Structure Database (ICSD).⁵² The supercell cell size for Na₃SbS₄ was set to 2 × 2 × 2 (for the cubic phase, this results to ∼14 Å cell edge dimension in three crystallographic directions). Optimization convergence was set to <1 meV per atom in energy and <0.01 eV Å⁻¹ in residual forces.

Thermodynamics calculation

Competing phases were referred from the Materials Project database via its application programming interface which also provided tools in this work to aid in phase diagram construction and analysis.^53,54 Na₃SbS₃, a compound which was separately reported in a previous work,³⁷ was included in the list of compound entries (resulting in a total of 61 phases) for the phase diagram construction.

The Gibbs free energy (G) at a given temperature (T) and pressure (P) was then calculated by:⁵⁵

G(T,P,N_Na,N_Sb,N_S) = H(T,P,N_Na,N_Sb,N_S) + PV(T,P,N_Na,N_Sb,N_S) − TS(T,P,N_Na,N_Sb,N_S), (1)

where H is the enthalpy, V is the volume, and S is the entropy. For condensed phases at 0 K, TS and PV terms were assumed to be negligibly small. Reference chemical potentials for the Na–Sb–S system were based on Na metal (μ⁰_Na), Sb metal (μ⁰_Sb), and S solid (μ⁰_S).

For (electro)chemical stability analysis, the grand potential phase diagram was formulated,⁵⁵ identifying the phase equilibria with respect to open Na species (μ_Na). The applied electrostatic potential ϕ is defined as:

μ_Na(ϕ) = μ⁰_Na − eϕ. (2)

In here, the decomposition energy of NSS at a given ϕ (ΔE_d(ϕ)) was determined using:^54–56

ΔE_d(ϕ) = E_equil(C_equil([x_Na,x_Sb,x_S],ϕ)) − E(Na₃SbS₄) − Δn_Na × (μ_Na(ϕ)), (3)

where E_equil(C_equil([x_Na,x_Sb,x_S],ϕ)) is the phase equilibria energy from the composition system C_equil([x_Na,x_Sb,x_S]) at ϕ, E(Na₃SbS₄) is the energy of NSS and Δn_Na is the number of Na exchanged between NSS and C_equil([x_Na,x_Sb,x_S]) during sodiation or desodiation.

Chemical stability analysis was also performed under “rich” and “poor” conditions of non-Na species in NSS, SE-related phases and possible cathode materials using the methodology in ref. 70. Chemical potentials for open Sb and open S species are referenced as well to μ⁰_Sb and μ⁰_S, respectively, which means that under Sb-rich (e.g., direct contact with Sb metal) and S-rich (e.g., direct contact with S solid) conditions, μ_Sb = 0 and μ_S = 0, respectively. Going towards Sb-poor or S-poor conditions means μ_Sb < 0 and μ_S < 0, respectively.

Elastic property calculation

Geometry optimizations were performed using stricter settings: kinetic energy cut-off of 700 eV, electronic self-consistency convergence of 10⁻⁷ eV, and at least 4 × 4 × 4 k-point. Optimized structures were used as initial cells to determine the elastic tensor. The symmetry-general least-squares extraction approach of elastic coefficients for strained materials was implemented, with elastic moduli including both contributions for distortions with rigid ions and from the ionic relaxations.⁵⁷ The elastic tensor was determined by carrying out 6 finite lattice distortions, after which the elastic constants were calculated from the resulting strain–stress relationship. For the bcc Na metal, experimental elastic data were used.⁵⁸

Defect calculation

The NSS unit cell was expanded into a 2 × 2 × 2 supercell (128 atoms). The defect formation energy (E_f,d) was calculated using the expression:where E_total,d is the total energy of the supercell with defects, E_{total, pristine} is the total energy of the pristine supercell, Δn_i is the number of atoms of element i that is added to or removed from the supercell (summed up for all N elements) to satisfy the charge balance condition, and μ_i is the chemical potential of species i.

Molecular dynamics calculation

For molecular dynamics (MD) calculations, a 2 × 2 × 2 c-NSS supercell with ∼14 Å cell edges was used. A cutoff kinetic energy of 400 eV was set and the k-point grid was fixed to 1 × 1 × 1. A revised PBE-GGA functional for densely packed solids (PBEsol) was used to describe electron exchange–correlation energy.⁵⁹

Different supercells were constructed for both undoped and halide-doped systems. In the doped models, 2 dopant atoms are introduced homogeneously at 2 S sites (see Fig. S1†) with charge compensation by Na vacancy (∼4%).

A comparison of the energetic stability between the tetragonal and the cubic NSS phase was performed and it showed the increasing stability of the latter with increasing Na vacancy concentration and halide doping content (see Fig. S2†). In addition, it was determined both by the present DFT analysis that the tetragonality of the NSS structure weakens with increasing number of Na vacancies and the structure becomes more cubic, this is consistent with an experimental study.⁶⁰ Specifically, there is a decrease in the ratio between the longer and the shorter lattice parameters of the tetragonal structure with Na vacancy introduction via halide doping (see Table S1†). Based on these results, the use of the cubic structure for the MD runs of halide-doped systems is reasonable and justified.

Standard pseudopotentials were employed for Na, Sb, S, F, Cl, Br, and I. The MD step size was set to 1 fs. Sampling temperatures were chosen in the range of 373–973 K, at every 100 K interval. To account for thermal expansion, an equilibration step was performed for 10 000 MD steps using a Langevin thermostat (NPT ensemble).⁶¹ The last 3000 MD steps were used to estimate the average lattice constant of the cubic cell at every trajectory-sampling temperature with 100 000 steps for the production run. The latter procedure used a thermostat in the Nose–Hoover formalism (NVT ensemble).⁶²

To characterize the sampled trajectories, the mean square displacement (MSD) of the atoms was calculated:

MSD = 〈[r⃑(t + τ) − r⃑(t)]²〉, (5)

where r⃑(t), t and τ are the atom position, time and lag time, respectively. Diffusion coefficient (D) calculation was based on the Einstein–Smoluchowski equation:⁶³where d is the dimensionality-related parameter related to the characteristic diffusion process. For a 3D ion conducting material, d = 3. The value of D can be estimated from the slope of the fitted linear function at the diffusive regime of the MSD profile.

Results and discussion

Crystal structure and electronic properties

The crystal structures of t-NSS and c-NSS are shown in Fig. 1. In t-NSS, isolated SbS₄ tetrahedral units (Wyckoff 2b Sb, 8e S) form a body-centered-type arrangement and 2 distinct Na sites (Na1: 4d, Na2: 2a) in distorted polyhedral units are positioned inside 6-S and 8-S cages. Two characteristic local pathway linkages form a contiguous 3D network for long range Na inter-diffusion: the ab-plane pathway with alternating Na1 and Na2 sites and the c-direction pathway which comprises only Na1 sites. In the case of c-NSS, only one Na Wyckoff site exists and the resulting pathway is 3D-symmetric.

Fig. 1 (a) Crystal structure of t-Na₃SbS₄ (P4̄2₁c) and (b) c-Na₃SbS₄ (I4̄3m). Highlighted are Na atoms and isolated SbS₄ tetrahedral units, linkage of Na Wyckoff sites in the ab-plane, linkage of Na Wyckoff sites in the c-direction, and 3D perspective-view of the Na pathway.

DFT-calculated lattice information for NSS is summarized in Table 1. Values are noted to be overestimated at most by 1.2% and 3.7% for the lattice parameters and volume, respectively, vs. experiment. Experimental bond distances are reproduced by DFT to within 2.3% (Table S2†). These discrepancies are error levels typical for the GGA-PBE functional.⁶⁴ From DFT total energy comparison, t-NSS was predicted to be more stable than the c-NSS, this is also consistent with experiment. Halide anion addition into the c-NSS structure results in a slightly expanded lattice cell by about 0.43–0.87% relative to its undoped form. Overall, the lattice constant increases with increasing halide anion radius (see Table S3†).

Table 1 Comparison between DFT results (in this work) and experiment (in parenthesis)³³ for the structural information on Na₃SbS₄

Phase	Lattice cell information				DFT total energy/eV per cell
	a, b, c/Å	Δ/%	V/Å^3	Δ/%
t-Na3SbS4 (P4̄21c)	7.23 (7.14), 7.23 (7.14), 7.37 (7.28)	1.2, 1.2, 1.2	385.3 (371.5)	3.7	−61.4313
c-Na3SbS4 (I4̄3m)	7.25 (7.17), 7.25 (7.17), 7.25 (7.17)	1.1, 1.1, 1.1	380.5 (368.6)	3.2	−61.3648

---

Based on the density of states (DOS) profile of t-NSS (Fig. S3†), the top of valence bands and the bottom of the conduction bands are found to comprise bonding and antibonding combinations of Sb and S, respectively. The bottom of the conduction band is split into: (i) S antibonding states and Sb 5s states for the lower band and (ii) S antibonding states and Sb 5p states for the upper band. The estimated GGA-PBE electronic band gap energy is 2.1 eV (we calculated at the HSE level a value of 3.1 eV), this is comparable with other sulfide-type SEs such as Li₁₀GeP₂S₁₂ (3.6 eV).⁶⁵ Usually, the band gap energy is used as a parameter for the SE electrochemical stability window (positive correlation) and electronic conductivity (negative correlation).

Thermodynamic phase stability, possible secondary phases

Using the formation energy data of competing phases in the Na–Sb–S system (see Table S4†), the convex hull results point to t-NSS as a ground-state phase (see Fig. S4†). Unless explicitly mentioned hereafter, thermodynamic-related analyses will be based on t-NSS. For (electro)chemical stability evaluation with cathode/anode contact and ion conduction property analysis, possible secondary/decomposition phases that may appear during material synthesis and battery cycling were also considered. This would enable the following: (i) the direct analysis of bulk-related properties and (ii) the qualitative analysis of inter-grain-related properties. During the material synthesis of NSS, chemical potentials of Na and S can be controlled (e.g., via Na₂S loss at high temperature to achieve ‘Na-poor’ synthesis conditions, or by explicit addition of excess Na₂S to realize ‘Na-rich’ conditions); this can lead to formation of other phases. Based on the Na–S grand potential phase diagram as shown in Fig. 2, possible secondary/decomposition phases were identified from the bounding phases of the NSS chemical potential stability region, these are Na₃SbS₃ (P2₁3), NaSbS₂ (C2/c), and Na₂S (Fm3̄m).

Fig. 2 DFT-calculated chemical potential (μ) range of thermodynamic stability for SE-related compounds evaluated in this study. The stability range is described by the upper bound μ₁ (in the rich-condition environment of a given chemical species) and lower bound μ₂ (in the poor-condition environment). Outside this μ range, a given compound will undergo decomposition. No stable phases exist in the shaded region. Note: μ_Na = 0, μ_Sb = 0 and μ_S = 0 are referenced to the chemical potentials of Na metal (μ⁰_Na), Sb metal (μ⁰_Sb), and S solid (μ⁰_S), respectively.

Cathode-side stability

Interface stability at the SE–cathode contact is an important issue that directly affects the cycling performance of solid-state batteries. It may be evaluated according to the thermodynamic tendencies of chemical species to migrate across the SE–cathode interface which could trigger the interface atomic rearrangement and/or interfacial decomposition, eventually leading to an increased interfacial resistance and poor capacity retention. Some reports have shown the formation of a reaction layer which may have been facilitated by species diffusion at the solid hetero-interface.^66,67 Other studies depict an interface region based on the space-charge layer effect in the scale of few-to-several atomic layers.⁶⁸ In here, the tendency for chemical species migration was qualitatively investigated using bulk chemical potential relationships to provide insight into the chemical driving forces that would cause the compositional change around the interface. For the cathode counterpart, the layered structures (Na[TM]O₂, where TM = {V, Cr, Mn, Fe, Co, Ni}) were considered here as contact materials with the electrolyte. To identify the possible reactions, the chemical potentials of the different species from the cathode and SE component were determined from element stoichiometry and bulk Gibbs free energy of relevant phases.

In the case of reductive and oxidative stability evaluation related to battery charge/discharge, it is noted that the local equilibrium of Na⁺, e⁻, and Na determines the electrode potential. It is assumed here that the local equilibrium may be established and the phase stability range of the SE compound can be defined (i.e., vs. Na⁺/Na).⁶⁹ This can be extended as well to other chemical species in the system. The phase stability range ([μ₁, μ₂]) can be described with respect to a specified open chemical species, where μ₁ and μ₂ denote bounds towards “rich” and “poor” conditions, respectively.^70–72 Outside [μ₁, μ₂], the equilibrium is pushed towards decomposition which forms the stable competing phases. When viewed at the cathode–SE contact, species migration is assumed to be likely if no common μ exists (for both the cathode and SE) that would allow for a possible equilibrium μ condition across the interface (i.e., checking if [μ₁,μ₂]_cathode∪[μ₁,μ₂]_SE = ∅). Otherwise, there would be no driving force for species migration (i.e., if [μ₁,μ₂]_cathode∪[μ₁,μ₂]_SE > ∅). For the cathode counterpart, experimentally confirmed compounds with layered structures were selected, such as NaCrO₂ (R3̄m).⁷³ Overall, 6 cathode compounds under pre-charged (i.e., full sodiation) and charged (i.e., full desodiation) conditions, as well as 4 SE/SE-related compounds (NSS, Na₃SbS₃, NaSbS₂, Na₂S) were investigated (total: 6 × 2 × 4 = 48 different interface–contact combinations).

Fig. 3a shows the [μ₁, μ₂] relationship for various cathodes and SE/SE-related compound pairings. The horizontal axis represents μ of the open TM species while the vertical axis represents the index/label for various layered Na[TM]O₂ cathode compounds. The horizontal lines for different compounds represent the [μ₁, μ₂] stability range. Additionally, for the different cathodes, the horizontal lines indicate the pre-charging condition (full Na content) while the corresponding horizontal bands are for the charging condition (full desodiation cases). An example for the assignment of μ₁ and μ₂ is shown in the top subplot for Ni species in NaNiO₂ (in red). Meanwhile, actual [μ₁, μ₂] values for open V, Cr, Mn, Fe, Co, and Ni species in NaVO₂, NaCrO₂, NaMnO₂, NaFeO₂, NaCoO₂, and NaNiO₂, respectively, are determined to be [0 eV, −1.79 eV], [0 eV, −4.01 eV], [−0.25 eV, −3.55 eV], [0 eV, −5.03 eV], [−0.08 eV, −2.78 eV], and [−0.36 eV, −2.15 eV], respectively. For NSS, μ ranges for open V, Cr, Mn, Fe, Co, and Ni species are [−3.27 eV, −∞], [−2.06 eV, −∞], [−1.62 eV, −∞], [−1.50 eV, −∞], [−1.01 eV, −∞], and [−0.84 eV, −∞], respectively. Based on the (non-)overlap in the μ ranges between the cathodes and SE/SE-related phases, the tendency for TM cation migration across the interface under battery pre-charging (or simple contact) and charging conditions is predicted to be unlikely, except for TM = V. The schematic summary is shown in Fig. 4a and detailed reactions are summarized in Tables S5 and S6.†

Fig. 3 (a) Chemical potential stability [μ₁,μ₂] range for TM species in Na[TM]O₂ (TM = {V, Cr, Mn, Fe, Co, Ni}) when in contact with t-Na₃SbS₄ solid electrolyte (SE) and other sulfide-type SE-related phases. (b) [μ₁,μ₂] range for X species (X = {O, Sb, S}) for various Na[TM]O₂–SE (SE-related) interface combinations. Horizontal bands are for same cathodes under charging condition (full desodiation).

Fig. 4 Schematic picture for (a) transition metal (TM) species and (b) S species tendency to migrate across the Na[TM]O₂–SE interface based on the chemical potential relationship as shown in Fig. 3a and b, respectively.

Fig. 3b shows the [μ₁, μ₂] relationship for open species X = {O, Sb, S}. For X = O, all 4 SE-related phases have common μ values with TM = {V, Cr, Mn, Fe} but not with TM = {Co, Ni} under the pre-charged cathode condition (i.e., full sodiation, horizontal lines in Fig. 3b). For X = Sb, overlap in [μ₁, μ₂] is noted for pairs NSS–TM = {V, Cr, Mn, Fe, Co}, SE = {Na₃SbS₃, NaSbS₂} – TM = {V, Cr, Mn, Fe }, and Na₂S–TM = {V, Cr, Mn, Fe, Co, Ni}. No common μ values appear for NSS–TM = Ni, and SE = {Na₃SbS₃, NaSbS₂} TM = {Co, Ni}. For X = S, species migration is predicted to be favorable between NSS and all Na[TM]O₂ compounds. Na₃SbS₃ has [μ₁, μ₂] overlap for TM = {V, Cr, Fe}, while for NaSbS₂ only overlaps for TM = {V, Cr}. Na₂S has overlap for TM = {V, Cr, Mn, Fe, Co}, but not for TM = Ni. These results are similar to the experimental observation for the Li-based cathode–SE interface, especially on S species migration tendency.^66,67 A schematic picture is shown in Fig. 4b for S species (see Fig. S5 in the ESI† for X = {O, Sb}). Under charging conditions (i.e., full desodiation, horizontal bands in Fig. 3b), the [μ₁, μ₂] ranges for X = O for the cathodes become largely narrow, their overlap with the [μ₁, μ₂] ranges of NSS/SE-related phases disappear (such as for TM = {V, Cr, Mn, Fe}), suggesting the increased chemical driving force for O species migration from the charged cathode to the SE side. The same trends are observed for X = {Sb, S}. These findings highlight the reactive nature of the charged cathode which can facilitate crossing of species and reactions at the cathode–SE contact. Detailed reactions are listed in Tables S7 and S8.†

The above results for the pre-charged cathode are found to be consistent with a previous DFT study. The result indicates the less favorable displacement reaction of S²⁻ by O²⁻ for NSS and pre-charged layered-type cathodes.⁷⁴ The contrast between NSS and the other electrolyte Na₃PS₄ was also confirmed, that is, the contact between Na₃PS₄ and layered-type cathodes highly favors PO₄ formation.⁷⁴ This difference between NSS and Na₃PS₄ may be further explained by the hard–soft acid–base theory, which emphasizes that Sb⁵⁺ (a soft acid) prefers to bind with S²⁻ (a soft base) while P⁵⁺ (a hard acid) prefers to bind with O²⁻ (a hard base).³²

Anode-side stability

The (electro)chemical stability of NSS is evaluated here by simulating conditions at the low (reduction, high μ_Na) and high (oxidation, low μ_Na) voltage regimes of a typical battery cycling process. The evaluation is assumed to be closely linked to Na insertion/de-insertion into/from the SE structure. Two possible scenarios could occur: (i) no reductive/oxidative decomposition or zero Na net loss for SE (i.e., Na⁺ ions are simply being shuttled between the anode and the cathode) or (ii) non-zero Na net loss for SE which is associated with irreversible side reactions (i.e., Na⁺ ions are trapped inside decomposition phases and become excluded from the shuttling process between the anode and the cathode). Table 2 summarizes the competing phases that would likely form during electrochemical cycling. Overall, narrow voltage windows were found for NSS, Na₃SbS₃, and NaSbS₂; calculated ranges are [1.60 V, 2.37 V], [1.18 V, 1.69 V], and [1.35 V, 2.43 V], respectively. These materials are predicted to react and decompose, both reductively (e.g., at <1 V) and oxidatively (e.g., at >2.5 V). Meanwhile, Na₂S shows a window of [0 eV, 1.94 V] which indicates stability at the low voltage regime. If it forms after the first few cycling steps or it already exists in the pre-cycled SE powder/pellet, it can form a passivating phase (electronic insulator) that may prevent further SE decomposition at the anode side. It is noted that the change in volume between NSS and its decomposition products (at the reductive regime) can be significant enough to possibly result in stress build-up and eventually, cracking and contact loss around the interface region.⁷⁴

Table 2 Summary of DFT-predicted reaction product phases related to the voltage stability window for Na with t-Na₃SbS₄ and related phases

Compound	Voltage/0 V vs. Na^+/Na	Reaction product phases
Na3SbS4	<1.60	Na3Sb, NaSb, Sb, Na3SbS3, Na2S
	>2.37	NaSbS2, Sb2S3, S
Na3SbS3	<1.18	Na3Sb, NaSb, Sb, Na2S
	>1.69	Na3SbS4, NaSbS2, Sb2S3, S
NaSbS2	<1.35	Na3Sb, NaSb, Sb, Na2S
	>2.43	Sb2S3, S
Na2S	>1.94	S

---

Anode-side dendrite formation tendency

Dendrite formation is another major hindrance towards developing batteries with metal anode design. In here, dendrite initiation/growth tendency on the Na metallic anode surface in contact with an SE material was evaluated using the Monroe–Newman 2D model.^75,76 This model evaluates the effects of general stresses on current density distributions at roughening interfaces. Starting from the typical reduction and deposition of Na⁺ ion as given by the reaction:

Na⁺ + e⁻ ⇌ Na, (7)

the time evolution of the Na metal surface growth can be tracked according to the equation:where x is the general normal direction towards the SE material, t is the time, e_z is the unit vector parallel to the general metal surface plane, e_n is the unit vector normal to the local curvature of the metal surface, i is the current density, V_Na is the molar volume of Na atom in the anode, z is the charge of Na⁺ ion, and F is Faraday's constant. The surface overpotentials η and i are related via the Butler–Volmer equation:where i₀ is the exchange current density, α_anode and α_cathode are the charge transfer coefficients for the anode and cathode, respectively, R is the gas constant, and T is the temperature. The infinitesimal change of electrochemical potential dμ for the electron at the deformed interface is coupled to pressure p by partial differentiation:

Taking the ratio between the deformed and ideal interface (using eqn (9)), the rate of change of the electron electrochemical potential (Δμ_e⁻) is derived as:^75,76

where i_deformed and i_ideal are current densities of the deformed and ideal interface, respectively. In terms of general stresses and molar volume of Na⁺/Na, Δμ_e⁻ is given by:where V_Na⁺ is the molar volume of Na⁺ ion in the SE material, γ is the surface tension of the Na metal anode, k is the mean curvature at the deformed interface, τ^Na_d and τ^SE_d are the deviatoric stresses in the Na metal anode and SE material, respectively, and Δp^Na and Δp^SE are the gage pressures in the Na metal anode and SE material, respectively. To evaluate electrodeposition stability, Δμ_e⁻ is analytically transformed into the Fourier space as:

Δμ_e⁻ = χRe{A e^ikx}, (13)

where χ is regarded as the electrodeposition stability parameter and A is the prefactor of the real part of the complex expression. This formulation has χ (in kJ mol⁻¹ nm⁻¹) expressed as a function of elastic properties, surface roughness, and molar volume quantities:where υ is the molar volume ratio of the Na⁺ ion in SE to the Na metal atom in the anode, G_Na is the shear modulus of the Na metal anode, G_SE is the shear modulus of SE, ν_Na is the Poisson's ratio of the Na metal anode, and v_SE is the Poisson's ratio of SE. The first, second, and third terms are contributions from surface tension, deviatoric stress, and hydrostatic stress, respectively, at the interface.

When χ > 0, electrodeposition becomes unstable (i.e., perturbation is in phase with current density), resulting in a faster deposition at the metal surface peaks rather than at the metal surface valleys (i.e., electron chemical potential is lower at the peaks than at the valleys) and thereby promoting Na dendrite formation (see Fig. 5a). When χ < 0, the reverse is true in which electrodeposition is predicted to be stable (i.e., perturbation is out of phase with current density), resulting in a faster deposition rate at the valleys than at the peaks of the deformed/rough Na metal surface (i.e., electron chemical potential is lower at the surface valleys than at the surface peaks), resulting in a more uniform deposition and no dendrite growth (see Fig. 5b). To determine χ, elastic data calculated by DFT and partial molar volumes of Na⁺ ion in the SE/SE-related compound structures from DFT-optimized structures are used.

Fig. 5 Schematic illustration for conditions leading to (a) Na dendrite growth and (b) stable Na electrodeposition, based on the relationship between electron chemical potential (μ_e⁻) and surface roughness (wavenumber k) at the Na metal surface (peak and valley).

Fig. 6a plots the contributions of the surface tension and general stresses to χ as a function of surface roughness k (in m⁻¹) for the Na–NSS interface. It indicates a stabilizing effect for surface tension (χ < 0) while a destabilizing effect for deviatoric and hydrostatic stresses (χ > 0); this is consistent with calculations performed for metal–SE interfaces for Li-based solid-state batteries.⁷⁶ Using the k value in a typical electrodeposition process (1 × 10⁸ m⁻¹),⁷⁵ the contributions of surface tension, deviatoric stress, and hydrostatic stress with NSS are −0.04 kJ mol⁻¹ nm⁻¹, +7.09 kJ mol⁻¹ nm⁻¹, and +15.68 kJ mol⁻¹ nm⁻¹, respectively, leading to a χ value of +22.73 kJ mol⁻¹ nm⁻¹ (see k-axis intercept in Fig. 6a inset). This suggests that Na dendrite formation is likely to occur during Na plating/stripping with metallic Na. No dendrite suppression is expected as well for other SE-related phases (Na₃SbS₃: +21.41 kJ mol⁻¹ nm⁻¹, NaSbS₂: +29.62 kJ mol⁻¹ nm⁻¹, Na₂S: +24.77 kJ mol⁻¹ nm⁻¹). A stability map using the ratio of molar volume quantities (V_m,Na⁺/V_m,Na) and shear moduli parameters (G_SE/G_Na) is displayed in Fig. 6b for various Na-based SE compounds, including non-sulfide-type ones (e.g., NaZr₂(PO₄)₃,²⁴ Na₃Zr₂Si₂PO₁₂,⁷⁷ Na₃AlF₆,⁷⁸ Na₃OBr²⁷). There are 3 regions in the plot: (i) density-driven stability region in the lower left, (ii) “region of instability” in the middle, and (iii) pressure-driven stability region in the upper right. All the considered SE compounds lie in region (ii) (i.e., V_m,Na⁺/V_m,Na < 0 and G_SE > G_Na). Overall, none of the investigated compounds are predicted to be stable against dendrite growth. Meanwhile, along the log(G_SE/G_Na) axis, NSS is noted to be closer to the frontier of region (i) than the non-sulfide types.

Fig. 6 (a) DFT-derived plot of electrodeposition stability parameter χ as a function of surface roughness wavenumber k. (b) Electrodeposition stability map plotted using the molar volume ratio of Na⁺ cation in the solid electrolyte (SE) and Na atom in the Na metal anode vs. their shear modulus ratio, respectively. Dashed curves in red and blue indicate frontiers for electrodeposition stability for a theoretically incompressible SE material (i.e., limiting Poisson's ratio, v_SE = 0.50) and an arbitrary material with v_SE = 0.21.

Intrinsic bulk defects

Table 3 lists the various defect complexes (in the Kroger–Vink notation) for the NSS compound. The lowest defect energies were noted for S Frenkel (1.26 eV), Na–S antisite (1.01 eV), and Na₂S Schottky (2.05 eV) (S Frenkel and Na–S antisite are visualized in Fig. 7). These values suggest that intrinsic disorder is maybe common in the material in its as-prepared form, in particular under the high-temperature synthesis condition.³⁰ Meanwhile, the Na-Frenkel model showed a defect energy of 0.15 eV, but this defect configuration naturally relaxes back to the normal Na sites, that is, Na and defect species redistribute according to the reaction:

Table 3 DFT-calculated intrinsic bulk defects in t-Na₃SbS₄

Reaction equation	E f,d/eV per formula unit
a Na Frenkel defect relaxes largely, all Na atoms occupy the Na Wyckoff sites back.
Na Frenkel (2nd nearest neighbor):	0.15^a
Sb Frenkel:	3.51
S Frenkel:	1.26
Na2S Schottky:	2.05
Na–Sb antisite:	3.33
Na–S antisite:	1.01
Sb–S antisite:	3.53

---

Fig. 7 Schematic illustration of (a) S Frenkel and (b) Na–S antisite from DFT-optimized supercells. Defect species such as , , , and and defect-related bond lengths are shown in red.

This is reflective of a Na pathway with energetically shallow local sites, favourable for fast Na diffusion. On the other hand, for S Frenkel and Na–S antisite, we found that after structural optimization the interstitial S actually prefers to form a chemical bonding with a corner S (S_corner) of a neighboring SbS₄ unit (see Fig. 7a and b) as evidenced by the shorter distance between S_int and S_corner (2.044 Å). Interatomic distance is also increased between S_corner and its center Sb (Sb_center) of the same SbS₄ unit (from 2.350 Å in the pristine model to ∼2.5 Å at post-relaxation), indicating Sb–S bond weakening due to S_int. The S_int–S_corner complex forms a basic bridging S unit which is similar to those observed in polysulfide compounds such as I4̄2d NaS₂ (in which the experimental S–S average bond length is 2.070 Å).⁷⁹ It is pointed out here that S–S bonding formation may facilitate the stabilization of Na vacancy defects via the charge compensation mechanism. Previous characterization by Rietveld refinement for the undoped NSS phase has indicated the presence of 5% vacancies in one of the Na Wyckoff sites.³²

Na ion transport property

For ground-state t-NSS, two distinct Na local pathways describe the 3D long-range inter-diffusion process (see Fig. 1a). Meanwhile, for metastable c-NSS which has a higher σ than t-NSS by a factor of ∼2, there is a 3D symmetric pathway of linked equivalent Wyckoff sites as shown by the DFT-MD 2D-projected trajectory in Fig. 8a. Similar to the strategy for other known metastable SE materials, c-NSS stabilization can also be realized via doping, for example by halide anions X (F, Cl, Br, I) at the S site. Similar to the intrinsic defect energies calculated for the undoped NSS phase, DFT doping energies for NSS by halide anions with introduction were also estimated to be sufficiently low (0.95–1.11 eV f.u.⁻¹) to be accessible during high-temperature synthesis.^30,80 Thus, we focus hereafter on c-NSS and investigate σ and effects of halide doping on it. For this, the Arrhenius relationship was employed:⁸¹where σ₀, T, E_a, and k_B are the conductivity prefactor, temperature, Na⁺ ion activation energy, and Boltzmann constant, respectively. To estimate σ and E_a, DFT-MD calculations were performed, with respect to dopant type and vacancy concentration. For this, 6 supercell models (2 × 2 × 2, ∼128 atoms) were constructed: 1-vacancy undoped (Na_2.9375SbS₄ as 1v-u model), 2-vacancy undoped (Na_2.875SbS₄ as 2v-u model), F doping (Na_2.875SbS_3.875F_0.125 as 2v-F model), Cl doping (Na_2.875SbS_3.875Cl_0.125 as 2v-Cl model), Br doping (Na_2.875SbS_3.875Br_0.125 as 2v-Br model), and I doping (Na_2.875SbS_3.875I_0.125 as 2v-I model).

Fig. 8 (a) 2D-projected atom trajectory density from 100 ps NVT-MD calculation (473 K) of the 1-vacancy c-NSS (I4̄3m) supercell model. (b) DFT-MD-derived Arrhenius plot for Na⁺ ionic conductivity at different temperatures of 1-vacancy undoped (1v-u), 2-vacancy undoped (2v-u), F-doped (2v-F), Cl-doped (2v-Cl), Br-doped (2v-Br), and I-doped (2v-I) supercell models. Experimental data for undoped Na₃SbS₄ are included for comparison.⁴¹

Fig. 8b shows the Arrhenius plot derived from temperature-dependent MSD profiles (see Fig S6†). We note that, for example, in the 1v-u model which is the closest representative of the reported undoped phase, the same order of magnitudes was calculated for E_a and room-temperature σ (σ_RT) vs. experimental data. Thus, comparison can be reasonably and systematically carried out based on our calculation results. Without doping, such as in 1v-u and 2v-u models, an increase from a single vacancy (∼2%) to two vacancies (∼4%) resulted in a decrease in the logσ vs. 1000T⁻¹ slope. The slope is directly correlated with E_a, the estimated values of 0.201 ± 0.042 eV to 0.159 ± 0.052 eV, respectively, reveal the importance of Na vacancies in activating the superionic conduction in NSS. Without vacancy, a negligible diffusivity was found as confirmed by a flat MSD trend (see Fig. S7†).

As previously pointed out, it is expected that E_a may be further modulated with respect to halide dopant type since the diffusion channel can be modified differently even at the same level of vacancy concentration, this is owing to the ionic size and electronegativity/ionicity difference among halide anions. To examine this, results for halide anion dopant radius r_X,⁸²E_a, and extrapolated σ_RT from DFT-MD runs were collected and summarized in Table 4 (experimental data included for comparison). The r_X parameter was initially chosen to check whether it has a direct correlation with E_a and σ_RT, noting that it may be correlated with the variation in the anion sublattice packing and thus, the diffusion pathway.

Table 4 Comparison between DFT-MD-based and experimental Na⁺ ion activation energy E_a and its room-temperature conductivity σ_RT

Composition	Dopant anion radius rX^a/Å	DFT-MD Ea/eV (expt.)	DFT σRT/10^−2 S cm^−1
a 6 nearest neighbors.^81 b Nominal Na3SbS4.^41
Na2.9375SbS4 (2% vacancy, 1v-u)	—	0.201 ± 0.042 (0.269)^b	0.828 (0.21)^b
Na2.875SbS4 (4% vacancy, 2v-u)	—	0.159 ± 0.052 (—)	0.220 (—)
Na2.875SbS3.875F0.125 (4% vacancies, 2v-F)	1.33	0.281 ± 0.078 (—)	0.178 (—)
Na2.875SbS3.875Cl0.125 (4% vacancies, 2v-Cl)	1.81	0.105 ± 0.017 (—)	4.489 (—)
Na2.875SbS3.875Br0.125 (4% vacancies, 2v-Br)	1.96	0.073 ± 0.007 (—)	10.240 (—)
Na2.875SbS3.875I0.125 (4% vacancies, 2v-I)	2.20	0.267 ± 0.029 (—)	0.253 (—)

---

DFT-MD E_a values reveal a concave up trend with increasing r_X as follows: 0.281 ± 0.078 eV with F doping, 0.105 ± 0.017 eV with Cl doping, 0.073 ± 0.007 eV with Br-doping, and 0.267 ± 0.027 eV with I doping. Correspondingly, the σ_RT trend shows a concave down behaviour: 0.178 × 10⁻² S cm⁻¹, 4.489 × 10⁻² S cm⁻¹, 10.240 × 10⁻² S cm⁻¹, and 0.253 × 10⁻² S cm⁻¹, respectively. Overall, Cl and Br doping appears to be promising as dopants for σ optimization. To study how r_X correlates with the changes in the diffusion channel, other physical parameters were also investigated such as Na site cage diameter (d₀) and site-to-site bottleneck diameter (d₁), as displayed in Fig. 9a. These parameters were determined by the Voronoi tessellation approach and using the DFT geometry-optimized structures.⁸⁴Fig. 9b and c display the ranges for d₀ and d₁ which are 2.1–2.7 Å and 1.60–1.67 Å, respectively. The undoped models (1v-u, 2v-u) show the smallest d₀ and largest d₁ values. For the F-doped model (2v-F), both its d₀ and E_a are the largest, but it has the smallest d₁ values relative to other models. It can be interpreted as the path bottleneck becoming smaller and thus E_a becoming larger with F doping. It is apparent that there is a critical value for d₁ (∼1.62 Å), below which E_a starts to increase significantly as shown in Fig. 9c (for 2v-F and 2v-I models). No monotonic decrease in E_a is determined with further increase beyond the noted critical value of d₁ (for 1v-u and 2v-u models). Therefore, it can be argued that, aside from inducing ionic conduction via Na vacancy creation, the dopant halide anion type modulates the diffusion channel for the moving Na⁺ ion by changing the Na site-to-site bottleneck size.

Fig. 9 (a) Schematic illustration of the Na diffusion channel isosurface in NSS. The isosurface is generated by the bond valence sum (BVS) approach.⁸³ Na sites labelled 1 and 2 are highlighted, with local pathway features described; d₀: Na site cage diameter and d₁: Na site-to-site bottleneck diameter. (b) DFT-E_avs. d₀ plot and (c) DFT-E_avs. d₁ plot for 1-vacancy undoped (1v-u), 2-vacancy undoped (2v-u), F-doped (2v-F), Cl-doped (2v-Cl), Br-doped (2v-Br), and I-doped (2v-I) supercell models.

Discussion

Considering that there are only a few practical SE material choices so far for ASS-NIB application, NSS may be expected as one of the base SE materials which satisfy the high-σ requirement for actual application. When it is incorporated in the actual solid-state cell, further improvement (e.g., via morphology control and interface engineering) is still necessary to satisfy the other metrics such as chemical/electrochemical compatibility vs. electrode materials and stability vs. dendrite initiation/growth during the anode-side Na plating/stripping process.

The tendency of S species in NSS to migrate across the interface with the cathode is predicted as one of the highly likely processes that would drive the chemical makeup variation and atomic rearrangement around the interface region. This predicted single-species-type migration may also form part as an elementary process of other experimentally suggested reactions such as mutual species exchange and defect-driven redistribution of chemical species across the cathode–SE contact.^66,67 Usually, these reactions are linked to resistance increase around the cathode–SE interface region.^85,86 One possible way to reduce this interfacial resistance is by introducing a thin oxide buffer layer which acts as a physical barrier that prevents interface reaction.^70,87,88

Similarly, the insertion of a physical barrier between the anode and SE component is also a possible approach to prevent or limit SE decomposition. For instance, NaCl was found to be able to slow SE degradation kinetics in NSS.³⁵ The formation of NaCl at the low-voltage regime (i.e., under high Na chemical potential) in Cl-doped NSS was verified to be energetically preferable in the present work according to DFT free-energy phase diagram calculations (see Table S9†). However, a trade-off with interface and grain boundary ion transport is to be expected, given the difference in σ of the inserted barrier compound. The use of an alloy metal anode may also address the aforementioned stability problem, since the chemical potential of Na in the alloy phase is relatively lower than in the pure Na metal phase. For example, instead of Na metal as the starting anode, Na₃Sb and NaSb alloys may be used instead which are more stable phases under low voltage conditions. However, a possible drawback of this approach would be volume expansion/contraction, particularly at high current densities. This volume change can cause poor cycle performance due to eventual contact loss. To address this issue on volume change, porous morphology structures or porous supports have been tested.⁸⁹ Meanwhile, to stop dendrite initiation at the Na metal anode surface, the electrodeposition stability map in Fig. 4b suggests two ways to achieve density-driven stability (i.e., region (i)): (i) make Na in the SE compound more densely packed than Na in the Na metal anode or (ii) appropriately pair the SE compound with an anode with lower effective molar volume for Na atom (e.g., alloy). In addition, intentionally increasing the Na metal surface roughness may also prevent dendrite formation, as shown in the total χ curve in Fig. 4 in which the value becomes negative at higher surface roughness k (e.g., at ∼6.41 × 10¹¹ m⁻¹). The density-driven stability (region (i)) may also be achieved via interface morphology control, by taking into account the anisotropy owing to the crystallographic orientation of Na metal and SE surfaces.⁹⁰ The stress–strain condition and electron density distribution around the anode–SE interface may also be controlled by modulating the cell stack pressure.⁹¹

Na ion transport in NSS is determined to be assisted by as revealed by DFT-MD results. The halide dopants modify the bulk σ value by: (i) controlling Na vacancy concentration and (ii) varying the size of the Na inter-site path bottleneck. An order of 10⁻² S cm⁻¹ is predicted to be accessible for Br-doped and Cl-doped NSS phases. A previous experimental report showed σ optimization by Cl doping in a material analogue Na_3−δPS_4−xCl_x (x = 0.2).⁹² Interestingly, the analysis indicated that this system is a case in which the Na vacancy is limited and Cl is abundant, rather than Na vacancy and Cl concentrations being equal (i.e., following a typical charge compensation for S²⁻ substitution with Cl⁻). It was suggested that, for charge compensation, the P oxidation state was changed from 5+ to 4+ (based on solid-state NMR analysis). It is noted here that 4+ is an uncommon electronic configuration for the P element, instead it is 5+ which is the natural oxidation state. This is in contrast with the present DFT work, in which the Na vacancy concentration depends on halide content for charge neutrality and there is no change in the oxidation state of Sb (such as from 5+ to 4+). So far, we found no experimental report yet that suggests the possible oxidation state change of Sb from 5+ to 4+ in NSS. Meanwhile, experiments have demonstrated a room-temperature conductivity maximum by Cl doping, the values are 9.0 × 10⁻⁴ S cm⁻¹ and 2.9 × 10⁻³ S cm⁻¹ for x_Cl = 0.05 and x_Cl = 0.0625, respectively.^60,93,94 These measurements are consistent with the present DFT-MD data for Cl doping, but the experimental doping concentrations are noted to be slightly lower. Although a higher conductivity is exhibited by our prediction, our simulated doping level (at x_Cl = 0.125) is apparently difficult to obtain experimentally because of decreased achievable SE pellet density.⁶⁰ The formation of new structure and/or secondary phases (with inherently lower conductivity) could also play a role in the observed decrease in experimental conductivity.⁶⁰

Out of all the bounding phases in the phase diagram for NSS, only Na₂S may offer a favorable grain boundary pathway between particles, as suggested by its relatively low Na⁺ ion migration energy (see Fig. S8 and S9†). Thus, suppression of other secondary phases is crucial to achieve high σ. In addition, intrinsic defects such as S Frenkel and Na–S antisite are expected to also affect σ since they may alter a long-range inter-diffusion process in the bulk, possibly even blocking the Na pathway (see Fig. 7). Control of these intrinsic defects should thus be considered as well, particularly during the actual material synthesis of NSS.

Conclusions

A comprehensive theoretical study based on DFT and advanced kinetic modelling was carried out to investigate the phase stability, (electro)chemical stability, electrodeposition stability, defect chemistry, and Na ion transport property of NSS SEs for ASS-NIB application. In the cathode side (e.g., using NaCrO₂), interface atomic rearrangement is highly likely to be facilitated by S chemical species diffusion. Meanwhile, NSS undergoes decomposition under the low voltage (Na insertion reaction) and high voltage (Na de-insertion reaction) conditions; the predicted electrochemical window is 1.60–2.37 V vs. Na⁺/Na. It was determined that Na dendrite formation is likely during the electrodeposition process on the Na metal anode surface with NSS as the contact SE material, highlighting the need for other strategies to be implemented such as control of morphology and cell stack pressure. Halide doping at the S site enables increasing the Na vacancy concentration and modulating the diffusion path bottleneck size to tune the conductivity and activation energy barrier. Among the dopants, Br⁻ and Cl⁻ anions are predicted to be the most promising ones, both with activation energies of about ∼0.1 eV at 4% Na vacancies. Secondary phases and intrinsic defects during NSS synthesis are possible hindrances for maximizing the conductivity. Na₂S, a possible secondary phase that could either exist in the as-prepared sample or could form during battery charge–discharge operation, may aid in the inter-grain transport because of its low Na⁺ ion migration energy.

Author contributions

R. J. and Y. T. conceived the study with the help of B. G., and H. K. T. Calculations were performed by R. J. Contributions were made by Y. T., B. G., and H. K. T. in relation to the analysis and discussion of thermodynamic reactions, chemical potential relationships, and ion transport. The manuscript was written by R. J., with input/comments/feedbacks from other authors.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

R. J. is thankful to the Japan Science and Technology Agency (JST) Precursory Research for Embryonic Science and Technology (PRESTO) program for the financial support. This research was also supported in part by JST through “Materials research by Information Integration” Initiative (MI²I)”, ALCA-SPRING Grant Number JPMJAL1301 and COI-NEXT Grant Number JPMJPF2016, by JSPS KAKENHI Grant Numbers JP19H05815 and JP21K14729, and by MEXT as Elements Strategy Initiative Grant Number JPMXP0112101003, Materials Processing Science project (“Materealize”) Grant Number JPMXP0219207397 and “Program for Promoting Researches on the Supercomputer Fugaku” Grant Number JPMXP1020200301. The calculations were performed on the supercomputers at NIMS (Numerical Materials Simulator), Nagoya University (CX400), Kyushu University (ITO) and the supercomputer Fugaku at the RIKEN through the HPCI System Research Project (project ID: hp210105). Visualization of crystal structures was made with the VESTA software.⁹⁵ The molecular dynamics trajectory in Fig. 8a was visualized using the VMD software.⁹⁶ Plots were generated using python Matplotlib.⁹⁷

Notes and references

1. C. Delmas, Adv. Energy Mater., 2018, 8, 1703137.
2. K. Kubota, J. Electrochem. Soc., 2015, 162, A2538–A2550.
3. S.-W. Kim, D.-H. Seo, X. Ma, G. Ceder and K. Kang, Adv. Energy Mater., 2010, 2, 710–721.
4. Y. Yuan, K. Amine, J. Lu and R. Shahbazian-Yassar, Nat. Commun., 2017, 8, 15806.
5. K. Mizushima, P. C. Jones, P. J. Wiseman and J. B. Goodenough, Mater. Res. Bull., 1980, 15, 783–789.
6. C. Delmas, C. Fouassier and P. Hagenmuller, Physica B+C, 1980, 99, 81–85.
7. T. A. Hewston and B. L. Chamberland, J. Phys. Chem. Solids, 1987, 48, 97–108.
8. A. K. Padhi, K. S. Nanjundaswamy and J. B. Goodenough, J. Electrochem. Soc., 1997, 144, 1188–1194.
9. S.-M. Oh, S.-T. Myung, J. Hassoun, B. Scrosati and Y.-K. Sun, Electrochem. Commun., 2012, 22, 149–152.
10. M. Ati, B. C. Melot, J.-N. Chotard, G. Rousse, M. Reynaud and J.-M. Tarascon, Electrochem. Commun., 2011, 13, 1280–1283.
11. R. Tripathi, T. N. Ramesh, B. L. Ellis and L. F. Nazar, Angew. Chem., Int. Ed., 2010, 49, 8738–8742.
12. S. Nishimura, M. Nakamura, R. Natsui and A. Yamada, J. Am. Chem. Soc., 2010, 132, 13596–13597.
13. P. Barpanda, G. Liu, C. D. Ling, M. Tamaru, M. Avdeev, S.-C. Chung, Y. Yamada and A. Yamada, Chem. Mater., 2013, 25, 3480–3487.
14. K. Ozawa, Solid State Ionics, 1994, 69, 212–221.
15. S. Komaba, W. Murata, T. Ishikawa, N. Yabuuchi, T. Ozeki, T. Nakayama, A. Ogata, K. Gotoh and K. Fujiwara, Adv. Funct. Mater., 2011, 21, 3859–3867.
16. W. J. H. Borghols, D. Lützenkirchen-Hecht, U. Haake, W. Chan, U. Lafont, E. M. Kelder, E. R. H. van Eck, A. P. M. Kentgens, F. M. Mulder and M. Wagemaker, J. Electrochem. Soc., 2010, 157, A582–A588.
17. H. Xiong, M. D. Slater, M. Balasubramanian, C. S. Johnson and T. Rajh, J. Phys. Chem. Lett., 2011, 2, 2560–2565.
18. E. Ferg, R. J. Gummow, A. de Kock and M. M. Thackeray, J. Electrochem. Soc., 1994, 141, L147–L150.
19. Y. Sun, L. Zhao, H. Pan, X. Lu, L. Gu, Y.-S. Hu and H. Li, Nat. Commun., 2013, 4, 1870.
20. Y. Wang, X. Yu, S. Xu, J. Bai, R. Xiao, Y.-S. Hu, H. Li, X.-Q. Yang, L. Chen and X. Huang, Nat. Commun., 2013, 4, 2365.
21. K. Xu, Chem. Rev., 2014, 114, 11503–11618.
22. H. Che, S. Chen, Y. Xie, H. Wang, K. Amine, X.-Z. Liao and Z.-F. Ma, Energy Environ. Sci., 2017, 10, 1075–1101.
23. F. G. Will and S. P. Mitoff, J. Electrochem. Soc., 1975, 122, 457–461.
24. A. K. Ivanov-Schitz and A. B. Bykov, Solid State Ionics, 1997, 100, 153–155.
25. Q. Ma, M. Guin, S. Naqash, C.-L. Tsai, F. Tietz and O. Guillon, Chem. Mater., 2016, 28, 4821–4828.
26. T. J. Udovic, M. Matsuo, A. Unemoto, N. Verdal, V. Stavila, A. V. Skripov, J. J. Rush, H. Takamura and S. Orimo, Chem. Commun., 2014, 50, 3750–3752.
27. H. Nguyen, S. Hy, E. Wu, Z. Deng, M. Samiee, T. Yersak, J. Luo, S. P. Ong and Y. S. Meng, J. Electrochem. Soc., 2016, 163, A2165–A2171.
28. Z. Zhang, E. Ramos, F. Lalère, A. Assoud, K. Kaup, P. Hartman and L. F. Nazar, Energy Environ. Sci., 2018, 11, 87–93.
29. S. S. Berbano, I. Seo, C. M. Bischoff, K. E. Schuller and S. W. Martin, J. Non-Cryst. Solids, 2012, 358, 93–98.
30. I.-H. Chu, C. S. Kompella, H. Nguyen, Z. Zhu, S. Hy, Z. Deng, Y. S. Meng and S. P. Ong, Sci. Rep., 2016, 6, 33733.
31. L. Zhang, K. Yang, J. Mi, L. Lu, L. Zhao, L. Wang, Y. Li and H. Zeng, Adv. Energy Mater., 2015, 5, 1501294.
32. H. Wang, Y. Chen, Z. D. Hood, G. Sahu, A. S. Pandian, J. K. Keum, K. An and C. Liang, Angew. Chem., Int. Ed., 2016, 55, 8551–8555.
33. A. Banerjee, K. H. Park, J. W. Heo, Y. J. Nam, C. K. Moon, S. M. Oh, S.-T. Hong and Y. S. Jung, Angew. Chem., Int. Ed., 2016, 55, 9634–9638.
34. L. Zhang, D. Zhang, K. Yang, X. Yan, L. Wang, J. Mi, B. Xu and Y. Li, Adv. Sci., 2016, 3, 1600089.
35. E. A. Wu, C. S. Kompella, Z. Zhu, J. Z. Lee, S. C. Lee, I.-H. Chu, H. Nguyen, S. P. Ong, A. Banerjee and Y. S. Meng, ACS Appl. Mater. Interfaces, 2018, 10, 10076–10086.
36. D. Zhang, X. Cao, D. Xu, N. Wang, C. Yu, W. Hu, X. Yan, J. Mi, B. Wen, L. Wang and L. Zhang, Electrochim. Acta, 2018, 259, 100–109.
37. L. E. Rush Jr, Z. D. Hood and N. A. W. Holzwarth, Phys. Rev. Mater., 2017, 1, 075405.
38. H. Wang, M. Yu, Y. Wang, Z. Feng, Y. Wang, X. Lü, J. Zhu, Y. Ren and C. Liang, J. Power Sources, 2018, 401, 111–116.
39. S. Xiong, Z. Liu, H. Rong, H. Wang, M. McDaniel and H. Chen, Sci. Rep., 2018, 8, 9146.
40. J. W. Heo, A. Banerjee, K. H. Park, Y. S. Jung and S.-T. Hong, Adv. Energy Mater., 2018, 8, 1702716.
41. A. Hayashi, N. Masuzawa, S. Yubuchi, F. Tsuji, C. Hotehama, A. Sakuda and M. Tatsumisago, Nat. Commun., 2019, 10, 5266.
42. T. Fuchs, S. P. Culver, P. Till and W. G. Zeier, ACS Energy Lett., 2020, 5, 146–151.
43. F. Tsuji, N. Masuzawa, A. Sakuda, M. Tatsumisago and A. Hayashi, ACS Appl. Energy Mater., 2020, 3, 11706–11712.
44. G. Kresse, Phys. Rev. B: Condens. Matter Mater. Phys., 1993, 47, 558–561.
45. G. Kresse and J. Furthmüller, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 54, 11169–11186.
46. J. P. Perdew and Y. Wang, Phys. Rev. B: Condens. Matter Mater. Phys., 1992, 45, 13244–13249.
47. G. Kresse and D. Joubert, Phys. Rev. B: Condens. Matter Mater. Phys., 1999, 59, 1758–1775.
48. P. E. Blöchl, Phys. Rev. B: Condens. Matter Mater. Phys., 1994, 50, 17953–17979.
49. J. P. Perdew, A. Ruzsinszky, G. I. Csonka, O. A. Vydrov, G. E. Scuseria, L. A. Constantin, X. Zhou and K. Burke, Phys. Rev. Lett., 2008, 100, 136406.
50. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865–3868.
51. H. J. Monkhorst and J. D. Pack, Phys. Rev. B: Solid State, 1976, 13, 5188–5192.
52. G. Bergerhoff and I. D. Brown, in Crystallographic Databases, F. H. Allenet al. (Hrsg.), International Union of Crystallography, Chester, 1987.
53. 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.
54. S. P. Ong, W. D. Richards, A. Jain, G. Hautier, M. Kocher, S. Cholia, D. Gunter, V. Chevrier, K. A. Persson and G. Ceder, Comput. Mater. Sci., 2013, 68, 314–319.
55. S. P. Ong, L. Wang, B. Kang and G. Ceder, Chem. Mater., 2008, 20, 1798–1807.
56. Y. L. Page, Phys. Rev. B: Condens. Matter Mater. Phys., 2002, 65, 104104.
57. G. Simmons and H. Wang, Single Crystal Elastic Constants and Cal-culated Aggregate Properties: A Handbook, Massachusetts Institute of Technology Press, Cambridge, 1971, p. 370.
58. C. Monroe and J. Newman, J. Electrochem. Soc., 2004, 151, A880–A886.
59. J. P. Perdew, A. Ruzsinszky, G. I. Csonka, O. A. Vydrov, G. E. Scuseria, L. A. Constantin, X. Zhou and K. Burke, Phys. Rev. Lett., 2008, 100, 136406.
60. F. Tsuji, S. Yubuchi, A. Sakuda, M. Tatsumisago and A. Hayashi, J. Ceram. Soc. Jpn., 2020, 128, 641–647.
61. W. G. Hoover, A. J. C. Ladd and B. Moran, Phys. Rev. Lett., 1982, 48, 1818–1820.
62. J. Nose, J. Chem. Phys., 1984, 81, 511–519.
63. P. Atkins and J. de Paula, Physical Chemistry, Oxford University Press, Oxford, 2006.
64. L. He, F. Liu, G. Hautier, M. J. T. Oliveira, M. A. L. Marques, F. D. Vila, J. J. Rehr, G.-M. Rignanese and A. Zhou, Phys. Rev. B: Condens. Matter Mater. Phys., 2014, 89, 064305.
65. Y. Mo, S. P. Ong and G. Ceder, Chem. Mater., 2012, 24, 15–17.
66. A. Sakuda, A. Hayashi and M. Tatsumisago, Chem. Mater., 2010, 22, 949–956.
67. R. Koerver, I. Aygün, T. Leichtweiß, C. Dietrich, W. Zhang, J. O. Binder, P. Hartmann, W. G. Zeier and J. Janek, Chem. Mater., 2017, 29, 5574–5582.
68. K. Takada, N. Ohta, L. Zhang, K. Fukuda, I. Sakaguchi, R. Ma, M. Osada and T. Sasaki, Solid State Ionics, 2008, 179, 1333–1337.
69. H. Yokokawa, Solid State Ionics, 2016, 285, 126–135.
70. W. D. Richards, L. J. Miara, Y. Wang, J. C. Kim and G. Ceder, Chem. Mater., 2016, 28, 266–273.
71. Y. Zhu, X. He and Y. Mo, ACS Appl. Mater. Interfaces, 2015, 7, 23685–23693.
72. Y. Zhu, X. Hea and Y. Mo, J. Mater. Chem. A, 2016, 4, 3253–3266.
73. S. Komaba, C. Takei, T. Nakayama, A. Ogata and N. Yabuuchi, Electrochem. Commun., 2010, 12, 355–358.
74. H. Tang, Z. Deng, Z. Lin, Z. Wang, I. H. Chu, C. Chen, Z. Zhu, C. Zheng and S. P. Ong, Chem. Mater., 2018, 30(1), 163–173.
75. C. Monroe and J. Newman, J. Electrochem. Soc., 2005, 152, A396–A404.
76. Z. Ahmad and V. Viswanathan, Phys. Rev. Lett., 2017, 119, 056003.
77. C. L. Tsai and H. Y.-P. Hong, Mater. Res. Bull., 1983, 18, 1399–1407.
78. M. Jansen, Angew. Chem., Int. Ed., 1991, 30, 1547–1558.
79. R. Tegman, Acta Crystallogr., Sect. B: Struct. Crystallogr. Cryst. Chem., 1973, 29, 1463–1469.
80. L. J. Miara, W. D. Richards, Y. E. Wang and G. Ceder, Chem. Mater., 2015, 27, 4040–4047.
81. R. J. D. Tilley, Defects in Solids, Wiley VCH, Hoboken, NJ, 10th edn, 2008.
82. R. D. Shannon, Acta Crystallogr., Sect. A: Cryst. Phys., Diffr., Theor. Gen. Crystallogr., 1976, 32, 751–767.
83. I. D. Brown, The Chemical Bond in Inorganic Chemistry - the Bond Valence Model, Oxford University Press, New York, 2002.
84. C. H. Rycroft, Chaos, 2009, 19, 041111.
85. B. Gao, R. Jalem, Y. Ma and Y. Tateyama, Chem. Mater., 2020, 32, 85–96.
86. H.-K. Tian, R. Jalem, B. Gao, Y. Yamamoto, S. Muto, M. Sakakura, Y. Iriyama and Y. Tateyama, ACS Appl. Mater. Interfaces, 2020, 12, 54752–54762.
87. N. Ohta, K. Takada, L. Zhang, R. Ma, M. Osada and T. Sasaki, Adv. Mater., 2006, 18, 2226–2229.
88. A. Sakuda, H. Kitaura, A. Hayashi, K. Tadanaga and M. Tatsumisago, J. Electrochem. Soc., 2009, 156, A27–A32.
89. L. Cao, T. Huang, Q. Zhang, M. Cui, J. Xu and R. Xiao, ACS Appl. Mater. Interfaces, 2020, 12, 57071–57078.
90. Z. Ahmad and V. Viswanathan, Phys. Rev. Mater., 2017, 1, 055403.
91. J.-M. Doux, H. Nguyen, D. H. S. Tan, A. Banerjee, X. Wang, E. A. Wu, C. Jo, H. Yang and Y. S. Meng, Adv. Energy Mater., 2020, 10, 1903253.
92. X. Feng, P.-H. Chien, Z. Zhu, I.-H. Chu, P. Wang, M. Immediato-Scuotto, H. Arabzadeh, S. P. Ong and Y.-Y. Hu, Adv. Funct. Mater., 2019, 29, 1807951.
93. H. Cao, M. Yu, L. Zhang, Z. Zhang, X. Yan, P. Li and C. Yu, J. Mater. Sci. Technol., 2021, 70, 168–175.
94. H. Gamo, N. H. H. Phuc, H. Muto and A. Matsuda, ACS Appl. Energy Mater., 2021, 4, 6125–6134.
95. K. Momma and F. Izumi, J. Appl. Crystallogr., 2008, 41, 653–658.
96. W. Humphrey, A. Dalke and K. Schulten, J. Mol. Graphics, 1996, 14, 33–38.
97. J. D. Hunter, Comput. Sci. Eng., 2007, 9, 90–95.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: 10.1039/d1ta07292g

---