Fast Na diffusion and anharmonic phonon dynamics in superionic Na₃PS₄

Abstract

The design of new solid electrolytes (SEs) hinges on identifying and tuning relevant descriptors. Phonons describe the atomic dynamics in crystalline materials and provide a basis to encode possible minimum energy pathways for ion migration, but anharmonic effects can be large in SEs. Identifying and controlling the pertinent phonon modes coupled most strongly with ionic conductivity, and assessing the role of anharmonicity, could therefore pave the way for discovering and designing new SEs via phonon engineering. Here, we investigate phonons in Na₃PS₄ and their coupling to fast Na diffusion, using a combination of neutron scattering, ab initio molecular dynamics (AIMD), and extended molecular dynamics based on machine-learned potentials. We identify that anharmonic soft-modes at the Brillouin zone boundary of the anharmonically stabilized cubic phase constitute key phonon modes that control the Na diffusion process in Na₃PS₄. We demonstrate how these strongly anharmonic phonon modes enable Na-ions to hop along the minimum energy pathways. Further, the quasi-elastic neutron scattering (QENS) measurements, supplemented with ab initio and machine-learned molecular dynamics simulations, probe the Na diffusivity and diffusion mechanism. These results offer detailed microscopic insights into the dynamic mechanism of fast Na diffusion and provide an avenue to search for further Na solid electrolytes.

---

Broader context

As part of the growing interest in renewable energy, the push for vehicle electrification motivates developing battery technologies based on safer and more abundant materials. While the dominant battery technology is currently Li-based, Li is a relatively rare element, present in large quantities in only few geographical locations, potentially posing a challenge for future growth in demand. Na constitutes a promising alternative as it is a low-cost and abundant chemical element worldwide. However, few materials exhibit sufficient Na-ion conductivity for battery applications. Hence, it is important to accelerate the discovery of new Na solid electrolytes (SEs), which requires a better fundamental microscopic understanding and the identification of descriptors relevant for fast Na-ion diffusion. Our investigation of atomic dynamics in Na₃PS₄, combining neutron scattering experiments and first-principles modeling, considers both the effects of Na vacancies and of the phonon modes of the host crystal. We identified anharmonic soft-modes that greatly facilitate the Na diffusion process in Na₃PS₄. Such anharmonic soft-modes could constitute a useful descriptor for designing new SEs. These results provide a microscopic understanding of how atomic vibrations facilitate the diffusion of Na ions through the material.

1 Introduction

The current push for safer rechargeable batteries with higher power densities is driving the search for solutions beyond the current Li-ion battery technology based on liquid or gel electrolytes, prompting investigations of all-solid designs leveraging solid-state electrolytes.^1–7 Besides safety concerns and performance limitations of liquid electrolytes, the high cost and limited availability of lithium are motivating research into batteries based on more abundant sodium.⁸ However, traditional Na-based solid electrolytes (SEs), such as β-alumina and NASICON compounds have achieved only limited ionic conductivity at room temperature. More recently, the thiophosphate Na ion conductor Na₃PS₄ and related compounds reached promising ionic conductivities near room temperature.^9–15 These Na superionic conductors (SICs) further offer promising tunability in their operating temperature by varying the host lattice composition. The atomistic mechanisms underlying fast ion conduction in solid electrolytes are currently the subject of intense scrutiny. In particular, a central question is how the vibrational properties of the crystal lattice, together with the presence of defects, may increase ionic diffusivity. The influence of the time-averaged structure and vacancies on the diffusion behavior has been studied in a range of SEs.^5,14–23 Studies of lattice dynamics in Na₃PS_4−xSe_x, Li-argyrodites, and Li-ionic conductors from LISICON and olivines families investigated the correlation between average phonon frequencies of mobile species and the hopping activation energy as well as the attempt frequency, and proposed that such an average phonon frequency could be a useful descriptor in the search for new SEs.^24–27 Yet, the characteristics of the host lattice that control Na or Li diffusion, in particular the influence of the flexibility of the lattice and of its large anharmonicity on the associated phonons, remain insufficiently understood to achieve predictive materials design.^27–30

Previous phonon studies in SEs drew useful correlations between lattice dynamics and diffusion, but remained largely limited to considerations of harmonic phonons, for example the average phonon frequency or density of states (DOS), thereby missing important considerations of mode-resolved dynamics and of the strong anharmonicity characterizing SICs.^31–33 In the harmonic picture of phonon theory, a Taylor expansion of the free energy surface is truncated at the second-order (harmonic), which results in independent phonon quasiparticles with well-defined energies.³⁴ However, this harmonic phonon picture is unable to account for properties such as diffusion, thermal expansion, thermal conductivity, or soft-mode phase transitions.³⁵ Higher-order anharmonic terms in the free energy need to be considered, and result in energy shifts and linewidth broadening of phonon spectra.³⁵ Because of their large-amplitude atomic motions and shallow potentials, superionic conductors tend to exhibit strong anharmonic effects.^31–33 Strong anharmonicity, as for example in the case of double-well potentials, may result in imaginary (unphysical) phonon frequencies within the harmonic approximation. Properly including higher-order anharmonic terms and finite temperature effects renormalizes phonon spectra, which can then be compared with experiments.^32,33,36,37

The relatively simple structure of Na₃PS₄ provides an opportunity to identify the structural and dynamical descriptors underlying its superionic behavior. The system crystallizes in a tetragonal phase (P4̄2₁c, α-Na₃PS₄) at low temperature and transforms to a superionic cubic phase (I4̄3m, β-Na₃PS₄) at 530 K^9,11 (Fig. 1).

Fig. 1 The crystal structure of (a) low temperature tetragonal phase (P4̄2₁c) and (b) high temperature superionic cubic (I4̄3m) phase of Na₃PS₄.

Further, Na₃PS₄ exhibits a high-temperature polymorph (γ-Na₃PS₄) for T ≥ 800 K, with an order-of-magnitude jump in ionic conductivity across the transition.¹⁹ The closely related compound Na₁₁Sn₂PS₁₂ (NSPS) also exhibits very high ionic conductivity.^38,39 Combining diffraction, impedance spectroscopy, and ab initio or classical molecular dynamics simulations, previous studies suggested that the underlying host dynamics in these materials play an important role in enabling the high Na mobility.^19,39,40 In γ-Na₃PS₄, the reorientational dynamics of PS₄ units (rotational jumps between symmetrically allowed PS₄ orientations) was proposed to be a key factor enabling high ionic conductivity.¹⁹ A classical MD study of Na₃PS₄ based on simple parameterized interatomic potentials also found a correlation between Na hopping and polyanion reorientations in the high-T γ phase, and a jump in Na-conductivity across the β–γ phase transition, but did not discuss phonons.⁴⁰ AIMD simulations of NSPS also suggested that the large Na conductivity is facilitated by PS₄ reorientation dynamics.³⁹ Ref. 19, 39 and 40 all invoked a “paddle wheel” mechanism (involving reorientation dynamics of structural subunits favoring diffusion). Yet, a fundamental understanding of the diffusion mechanism in terms of the mode-resolved phonons of the parent structure (wave-vector and branch-resolved eigenvectors), and of the importance of lattice potential anharmonicity, remains elusive.

Both NSPS and γ-Na₃PS₄ possess an open framework structure with a large volume per atom, in comparison to conventional solid ionic conductors.^19,39 However, few materials (such as Na₁₁Sb₂PS₁₂ and Li₆PS₅I) with larger volumes still exhibit limited diffusion coefficients compared to other family members.^39,41,42 Thus, an open-structure with large volume could be a necessary but non-sufficient condition to facilitate fast diffusion. Apart from the structure of the host lattice, its dynamics may also contribute to cation diffusion. In particular, the presence of low-energy modes (combining wiggling and librational motions of structural subunits) and their coupling with the diffusing ions may facilitate the diffusion process.^{19,39,40,43,44} Such low-energy phonons could be an important descriptor of ionic conductivity and need to be investigated thoroughly. Further, the traditional quasiharmonic phonon picture of conventional crystalline materials breaks down in superionic materials, owing to very strong anharmonic effects and dynamic disorder.^31–33 Thus, temperature-dependent measurements and simulations investigating phonon anharmonicity in superionic conductors are critical, yet remain scarce.^24,25,27

In this article, we report our investigations of atomic dynamics in Na₃PS₄ across a wide temperature range, considering on an equal footing the effect of Na vacancies and the phonons of the overall crystal structure, and their coupling, on the diffusive dynamics of Na ions. We use a combination of inelastic neutron scattering (INS) and quasi-elastic neutron scattering (QENS) experiments to obtain a full picture of the atomic dynamics. QENS probes the time and length scales characteristic of the diffusion process, while INS provides access to the phonon dynamics of the crystal lattice over a wide energy range and across the entire Brillouin zone (BZ). We performed temperature-dependent QENS and INS measurements of Na₃PS₄ across the phase transitions (100 K to 600 K) to probe the Na diffusion and the role of low-energy (E) phonons in the fast diffusion process. To rationalize our measurements and gain atomistic insights, we further performed ab initio molecular dynamics (AIMD) simulations, and extended these simulations to long time scales using machine-learned neural-network potentials. AIMD simulations are known for their high accuracy in predicting the thermodynamics properties but remain constrained by their high computational cost. Using machine-learning, we trained the surrogate potential against AIMD trajectories, which was then used to extend MD simulations to length and time scales several orders of magnitude larger than possible with AIMD while retaining high accuracy (“machine learned molecular dynamics”, MLMD).

Our experiments and simulations show that strongly anharmonic low-E phonon modes are already present in the normal phase and facilitate the Na-ion diffusion in the superionic phase. Specifically, facile Na jumps are enabled by the large amplitudes of low-E floppy dynamics of PS₄ subunits coupled with Na translations in 1D channels. We trace the key dynamic mechanism to zone-boundary soft modes (anharmonic low-energy modes) in the anharmonically stabilized cubic phase, which enable the fast ion conduction. These insights provide a potential new type of dynamic descriptor in screening for future solid electrolyte candidates, and highlight the importance of anharmonic effects in the lattice dynamics.

2 Materials and methodology

2.1 Neutron scattering

We performed the INS measurements using the high-flux ARCS time-of-flight⁴⁵ spectrometer at the Spallation Neutron Source (SNS) at Oak Ridge National Laboratory (ORNL). A polycrystalline sample of Na₃PS₄ (mass ∼ 6 grams) was encased in a 3/8′′ cylindrical aluminum sample holder and mounted in a high-temperature cryostat.

The INS data were collected using three incident neutron energies, E_i = 20, 40, and 100 meV. The measurements with E_i = 20 meV provide high resolution to resolve low-E peaks and accurately determine their energy shift with temperature, while the higher E_i data enable a survey covering the whole phonon spectrum, extending up to about 70 meV. The background from the empty aluminum holder was measured with the same settings at each temperature and subtracted from the measurement of the sample. The Mantid software was used for data reduction.⁴⁶

The dynamical structure factor, S(Q,E), contains comprehensive information on atomic dynamics. The neutron-weighted phonon DOS (g⁽ⁿ⁾(E)) in the incoherent approximation was obtained from measured S(Q,E) as described in ref. 47, after corrections for multiphonon and multiple scattering:

where Q and E is the momentum and energy transfer between neutron and sample, while + or − sign corresponds to the energy loss or gain of the scattered neutrons, respectively. T, k_B and A represent temperature, Boltzmann's constant, and a normalization constant respectively. The quantity between 〈⋯〉 represents a suitable average over Q values at a given energy, and 2W(Q) is the Debye–Waller factor averaged over all atoms.

The calculated g⁽ⁿ⁾(E) was obtained from the partial DOS (g_j(E)) as:

The weighting factors 4πb_j²/m_j (in units of barns per amu) are: 0.143, 0.107, and 0.032 for Na, P, and S, respectively, using neutron scattering lengths from ref. 48. The partial DOS, g_j(E) were obtained from MLMD simulation using eqn (6) described in the next section.

The QENS measurements were performed on the same powder sample as the INS experiment, using the BASIS backscattering spectrometer⁴⁹ at the SNS, ORNL. We used Si(111) analyzers with neutron wavelength of 6.4 Å and a chopper frequency of 60 Hz. In this configuration, the spectrometer provides an energy transfer range of −100 < E < 100 μeV with an energy resolution of 3.4 μeV, and momentum transfer coverage of 0.2 < Q < 2.0 Å⁻¹. The data were normalized and corrected for detector efficiency using a vanadium standard. A top-loading closed-cycle refrigerator with hot stage was used to measure at temperature T = 100 K, 300 K, 400 K, 500 K, and 600 K. The 100 K data do not show any broadening beyond the instrumental intrinsic width, as the Na dynamics are slow enough to be considered frozen, which were thus used as a resolution function for the analysis of QENS spectra at higher temperatures. The data from various detectors were grouped into Q bins of width 0.2 Å⁻¹. The data reduction and analysis were performed using the software packages Mantid⁴⁶ and Dave.⁵⁰ We fit the S(Q,E) from QENS using the following model:

where δ(E) in the square bracket represents the elastic signal, and the second term is the Lorentzian corresponding to quasi-elastic contribution with a half-width at half maximum Γ(Q), and R(Q,E) is the resolution function, for which we used the 100 K data. For each Q bin, the background function C(Q,E) is chosen to be a linear function in energy to reproduce the flat background from the sample environment and instrument. A and B are adjustable amplitudes of the elastic and quasi-elastic signals.

The Q dependence of Γ was analyzed using the Chudley–Elliot (CE) model for jump-diffusion:³¹

where d and τ are the average jump-length and residence time, respectively. Further, in the low-Q limit, the diffusion constant D is given by:

D = d²/6τ (5)

2.2 Computational modeling

Lattice dynamics. We performed first-principles lattice dynamics and AIMD simulations within the framework of plane-wave density functional theory (DFT), as implemented in the Vienna ab initio simulation package (VASP).⁵¹ The calculations used the projected augmented wave (PAW)⁵² with generalized gradient approximation (GGA) parametrized under the Perdew, Becke, and Ernzerhof (PBE) scheme.⁵³ Harmonic phonon calculations were performed in both the tetragonal and cubic phases within the small displacement approach implemented in Phonopy.⁵⁴ We used 2 × 2 × 2 supercells (128 atoms) and displacements amplitudes of 0.01 Å. The total energy and atomic forces were calculated on 14 (4) distinct displacement configurations in the tetragonal (cubic) phase. We used a plane-wave kinetic energy cutoff of 400 eV and on a 4 × 4 × 4 k-point grid generated by the Monkhorst–Pack⁵⁵ method. The self-consistent convergence threshold for electronic minimization was set to 10⁻⁸ eV. The anharmonic phonon dispersions at 600 K were calculated using renormalized force-constants derived from AIMD forces at 600 K using the software ALAMODE.⁵⁶ A second-order force constant cutoff of 7 Å was used.

Molecular dynamics. AIMD simulations were performed using a 2 × 2 × 2 supercell of the conventional cubic cell (128 atoms). Single k-points at Γ point and an energy convergence of 10⁻⁶ eV were used. We employed the NVT ensemble and ran the AIMD simulation for 12 ps with a time step of 2 fs. The temperature of the system was controlled by a Nosé–Hoover thermostat with a time constant of 0.1 ps. We equilibrated the system for 2 ps and the next 10 ps data for the production run. The DEEP-MD code⁵⁷ was used to train a surrogate neural-network force-field based on machine learning. Subsequently, long MLMD simulations were performed with this surrogate force-field using LAMMPS.⁵⁸ The surrogate force-field was benchmarked against AIMD results (pair distribution function, phonon DOS, and mean squared displacement (MSD)), showing excellent agreement (Fig. S9, ESI†). The MLMD simulations were performed on a 3 × 3 × 3 supercell (432 atoms) in pure and Na-vacant configurations with NVT (NVT-MLMD) and NPT (NPT-MLMD) ensembles. Temperature and pressure were controlled by a Nosé–Hoover thermostat and barostat with time constants of 0.1 ps and 1 ps, respectively. The trajectories were computed up to ∼1000 ps (using 1 fs time steps). These simulations provide momentum and energy resolutions of 0.2 Å⁻¹ and 1 μeV, respectively, essential to access the QENS spectrum through simulations. Further, to check the role of different host dynamics (e.g. rotation and wiggling of PS₄) on Na diffusion, we performed additional MLMD simulations with PS₄ units as rigid bodies. We investigated the effect of wiggling dynamics on Na diffusion by restricting the rotational degree of freedom of PS₄ tetrahedral units. In constrained MLMD simulations, we used a smaller supercell of dimension 2 × 2 × 2, due to the restriction of the maximum allowed numbers of rigid bodies (PS₄) in LAMMPS.⁵⁸ The choice of Na-vacancy concentrations (0.6% and 2%) is based on the minimum allowed Na vacancy concentration in the respective supercell dimensions. Our constrained-MLMD and NEB simulations on a 2 × 2 × 2 supercell allow a minimum of 2% Na vacancy concentration, while QENS and other simulations on a 3 × 3 × 3 supercell were performed with 0.6% and higher Na vacancy concentrations.

The partial (g_j(E)) and total phonon DOS (g(E)) were obtained from MD simulations via the Fourier transform of the velocity auto-correlation function as:

where v_j(t) is the instantaneous atomic velocity of element j at time t and 〈…〉 indicates an ensemble average over atoms.

The MSD(〈u_j²(t)〉) of element j at time t is given by:

where N_j is the total number of atoms of element j in the simulation cell, r_j(t) is the position of jth atom at time t and 〈…〉 represent the ensemble average. The isotropic diffusion coefficient of jth species at temperature T was estimated using the Einstein relation, i.e., from the time dependence of MSD over a period of τ:

We also computed the angular autocorrelation function, ς(t) of PS₄ tetrahedra, defined as:⁴⁰

ς(t′) = 〈r̂(t)·r̂(t + t′)〉 (10)

where r̂(t) is a unit vector connecting P and S within a tetrahedron at time t′, and 〈…〉 representing an average over multiple time origins.

3 Results and discussion

3.1 Phonon DOS and low-energy modes

Fig. 2(b) shows the resulting g⁽ⁿ⁾(E) and its temperature evolution, focusing on the low-E range. The DOS measured over a wider E range is shown in Fig. S1 (ESI†). As can be seen in Fig. 2(b), g⁽ⁿ⁾(E) exhibits peaks at very low energies E ∼ 6 meV at 100 K, which is surprising considering the low mass of the atoms in the compound. From our MLMD simulations, we extracted the species-decomposed partial DOS (Fig. S2, ESI†). We can see that the Na dynamics dominate g⁽ⁿ⁾(E) contributions below 30 meV, while the P and S contribute to the entire spectral range. It is unusual for a very light atomic species to exhibit such low-frequency vibrations (ω = (K/M)^1/2) with M the mass and K a characteristic force-constant), revealing the extremely soft bonding of the Na ions. The low-E peaks are also seen to soften (shift to lower energy) and broaden drastically on warming, revealing the strong anharmonicity of the interatomic potential, both in INS and MD simulations (Fig. 2(b and c)). We also observe some softening and more moderate broadening of other high energy peaks in the DOS on warming, indicating an overall softness and anharmonicity of the structure. The softening and broadening of low-E modes on warming is even clearer in the measured S(E) (Fig. 2(a)). We can also see that the low-E modes in Na₃PS₄ become strongly damped at high T, extending towards the elastic line in a quasielastic fashion (Lorentzian peak centered at E = 0). The large softening and broadening of these modes reflects that Na vibrations occur in strongly anharmonic potential wells and that the vibrational dynamics are strongly perturbed and overdamped when the rate of Na diffusion hops increases at higher temperatures, leading to a breakdown of phonon quasiparticles for Na-dominated modes, as detailed below.

Fig. 2 (a) Temperature-dependent INS measurements of the Q-integrated (1.0 Å⁻¹ ≤ Q ≤ 4.0 Å⁻¹) dynamical structure factor (S(E)) of Na₃PS₄, measuring the ARCS spectrometer with E_i = 20 meV. (b) Temperature-dependent neutron-weighted DOS derived from the INS measurement of S(Q,E). (c) Simulated neutron-weighted DOS obtained from the NVT-MLMD simulations at different temperatures. The measured and simulated DOS are normalized to unity over the energy range 0–15 meV.

Recently, Famprikris et al. reported temperature-dependent inelastic neutron scattering (INS) spectra, Raman spectra, and diffraction measurements above room temperature in Na₃PS₄. Their study focused on the nature of phase transitions, the local structure in different Na₃PS₄ polymorphs, and thermal expansion.²³ They suggested that low-energy modes below E ∼ 25 meV contribute to Na-hopping, but did not present spectra resolved to very low-energies. While our INS spectra agree well with INS measurements reported by Famprikris et al., our measurements with lower incident neutron energy (E_i = 20 meV) enabled us to critically resolve the key low-energy phonons below ∼7 meV that evolve into the quasielastic signal at high-T. Our QENS measurements of very low energy dynamics (E < 0.5 meV) further probe the Na hopping process. Combined with our anharmonic lattice dynamics simulations, the INS and QENS reveal the coupling of Na hoping with specific soft anharmonic phonons of the parent structure.

3.2 Diffusion pathway and phonon projection

To further investigate the role of phonons in Na diffusion, we performed nudged elastic band (NEB) simulations and projected these onto phonon eigenvectors. We start by considering the phonon dispersions in both cubic and tetragonal phases of Na₃PS₄ computed within the quasi-harmonic approximation (QHA), see Fig. 3(a) and Fig. S5 (ESI†). We observe stable phonon modes in the tetragonal phase, while the cubic phase exhibits unstable phonons (imaginary frequencies represented as negative values) revealing that it is dynamically unstable at T = 0 K in the harmonic approximation.

Fig. 3 (a) Calculated phonon dispersions at T = 0 K under QHA (dotted blue lines) and at 600 K using anharmonically renormalized force-constants (solid red lines) in the cubic phase of Na₃PS₄. (b) The eigenvector of the lowest unstable phonon mode at H-point (1/2 1/2 1/2) in cubic-Na₃PS₄. (c) Calculated 〈u²〉 as a function of phonon energy E averaged over entire Brillouin zone at 530 K in the cubic phase of Na₃PS₄. The modes below 10 meV significantly contributed to Na vibrations and PS₄ librations. (d) The calculated energy profile of the lowest unstable phonon mode at H-point (1/2 1/2 1/2) in cubic-Na₃PS₄. The vertical dashed line represent the configuration where all the Na atoms are at 6b sites (e) NEB calculation of energy barrier profile for Na migration between 6b crystallographic sites in stoichiometric phase (solid blue line) and 2% Na-vacant (dotted orange lines). The lowest phonon energy profile near 6b shown in (d) closely follows the stoichiometric NEB profile near 6b. (f) The calculated Na phonon DOS from NVT-MLMD in unconstrained, frozen-host (with only Na atomic degree of freedom) and frozen-rotation (PS₄ rotational degree are frozen, while PS₄ translations are allowed) simulation environment at 600 K. Constraining PS₄ dynamics results in suppression of low-E modes around ∼3–5 meV.

The maximum instability is observed at the H point, at (1/2, 1/2, 1/2) in reciprocal lattice units (rlu). The eigenvector of this mode, shown in Fig. 3(b), clearly exhibits librations of PS₄ units combined with Na translation along 〈100〉 channels. From an amplitude mode analysis,⁵⁹ we find that this H-point instability is the primary order parameter for the cubic–tetragonal phase transition. The unstable mode corresponds to a double-well potential for the potential energy along this eigenvector, characteristic of soft-mode instabilities (Fig. 3(d)). This indicates that symmetric PS₄ libration dynamics in the cubic phase freeze in with a finite amplitude in the tetragonal phase below T_c. The energy required to unlock these libration dynamics corresponds to the energy barrier between minima of the double-well (50 meV ≃ 580 K), which is fairly close to the observed transition temperature in Na₃PS₄.¹¹ On warming above T_c, anharmonic renormalization stabilizes phonons in the cubic phase, as shown in our finite-T phonon dispersion simulations (solid red lines in (Fig. 3(a)). The kinetic energy gain above T_c enables PS₄ units to oscillate between double-well minima and stabilizes the symmetric H-point mode, but its energy remains low, just a few meV.

At high T, such low-energy modes lead to large vibrational amplitudes, particularly Na translations and PS₄ librations and wiggling. This is shown by the decomposed MSD of individual elements as a function of phonon energy in Fig. 3(c). The MSD of atom j contributed by phonon of energy E was obtained as:

where g_j(E) and m_j are the partial DOS and mass of atom j. As can be seen in Fig. 3(c), the modes below 7 meV are dominated by Na dynamics. We also observed significantly larger MSD for S than for P, confirming the librational dynamics of PS₄ tetrahedra at low-E, although the P MSD also indicate wiggling motions of PS₄ units. At elevated T, the large population of these modes enhances the Na hopping probability from site to site along the [100] channels.

Importantly, the motions of PS₄ units widen the gateway for Na hopping, as we now detail. We have identified the minimum (near 12d site) and maximum diameters (near 6b) along the Na diffusion channels and their dependence on PS₄ motions (Table 1). We find that PS₄ motions indeed increase the channel diameter thus facilitating Na migration. This effect is similar to prior findings in the related compound NSPS.³⁸ Hence, cooperative effects between Na translation and PS₄ dynamics, captured in the anharmonic low-E phonon modes at the H-point, are critical for the fast Na hopping.

Table 1 The calculated minimum (R_min) and maximum (R_max) channel diameter and channel volume for Na-diffusion in Na₃PS₄ at equilibrium and saddle point configurations of NEB images⁶⁰

	(Rmin)/(Rmax) (Å)	Channel volume (Å^3)
At equilibrium	1.63/2.13	146.2
At saddle point	1.68/2.22	146.6

---

Our MD simulations further confirm that Na diffuses via jumps between 6b Wyckoff sites along the interconnected 1D channels. Hence, the most probable minimum energy pathway (MEP) is across 6b sites via the intermediate 12d sites (halfway between 6b sites along [100]). NEB calculations provide an estimated E_b ∼ 0.3 eV. In addition, we also find that the presence of Na-vacancies (∼2.0%) significantly lowers E_b, to ∼0.2 eV. By decomposing the NEB path onto the phonon eigenvector basis (ζ_NEB), we find that for small distortion amplitudes around 6b sites, the NEB profile and ζ^6b_NEB closely follow the H-point phonon potential and eigenvector (ζ^H_phonon) (Fig. 3(d and e) and Table S1, ESI†). However, the NEB path near the 12d site (ζ^12d_NEB) deviates significantly from ζ^H_phonon, indicating either that: (i) at higher distortion amplitude ζ^H_phonon changes significantly, and/or (ii) ζ^12d_NEB (Na pathways) cannot be represented by ζ^H_phonon. To understand this further, we performed calculations including the large distortion amplitude corresponding to the saddle-point configuration (moving 3 Na atoms from 6b site to 12d in a conventional unit cell). Strikingly, the phonon eigenvector with saddle-point configuration follows the ζ^12d_NEB (Table S1, ESI†). This further confirms that these low-E phonons indeed constitute precursors of the Na-hopping motions.

3.3 Hopping dynamics and diffusion mechanism

The temperature evolution of experimental QENS spectra at Q = 1.4 Å⁻¹ is shown in Fig. 4(a). A characteristic broadening is seen upon warming, which can be attributed to Na diffusion, as supported by our simulations. In particular, we find a negligible contribution of PS₄ dynamics to the QENS signal, see below. The Q dependence of QENS width contains rich information about the hopping process. The fits of QENS spectra are shown in Fig. S6 (ESI†). The fitting was performed over a Q range of 0.2 to 1.5 Å⁻¹. The results for Γ vs. Q at 600 K is shown in Fig. 4(b). The Q dependence of Γ follows the CE model well, and the estimated jump length is 3.6 Å. The estimated diffusion coefficient D is on the order of ∼10⁻⁶ cm² s⁻¹ at 600 K, typical of superionic conductors (Table 2).

Fig. 4 (a) Temperature dependent QENS spectra of Na₃PS₄, measured at Q = 1.4 Å⁻¹ (integrated over δQ = 0.2 Å⁻¹). (b) The measured Q-dependence of Γ (orange markers) fitted with Chudley–Elliot jump-diffusion model (solid line). (c) The calculated Γ obtained by fitting the Lorentzian to simulated S^Na_inc(Q,E) (dots) using NVT-MLMD trajectories with 0.6% Na-vacancy in Na₃PS₄ and Chudley–Elliot fit (solid line).

Table 2 The Na diffusion characteristics (d, τ and D) in Na₃PS₄ estimated from a Chudley–Elliot jump-diffusion model using the QENS measurement and simulated S(Q,E) (with 0.6% Na-vacancy and NVT ensemble) at 600 K

	d (Å)	τ (ps)	D (×10^−6 cm^2 s^−1)
QENS	3.6	246	0.9
MLMD	3.6	237	0.9

---

We performed the AIMD simulations at 100 K to 1000 K using an NVT ensemble. The calculated MSD for the stoichiometric compositions Na₃PS₄ did not show any signs of Na diffusion, as seen in Fig. S8 (ESI†), in agreement with previous investigations.¹⁴ Upon introducing Na vacancies in simulations, fast diffusion occurs in Na_3−xPS₄ (Fig. S8, ESI†). This result is intuitive, since fully filled 1D channels are expected to block the diffusion path. Hence, some vacant Na sites are needed to enable diffusion. To analyze the QENS spectra, we employed MLMD to reach sufficient long simulation times (see details in Materials and methodology section). This enables us to achieve the required energy and Q resolution to simulate the QENS spectra, which was impractical with AIMD simulations. We also note that AIMD trajectories of limited duration can result in significant errors in estimates of diffusion coefficients due to limited statistics of jump events, especially at low T.⁶¹

Further, we use MLMD simulations to assess the contributions to the QENS signal and resolve the element-specific coherent and incoherent contributions. While separating these components based on measurements alone is not practical. We identified that mobile Na atoms mainly contribute to the QENS spectra. The measured intensity in neutron scattering experiments is the sum of coherent and incoherent components:⁶²

where bⁱ_coh and bⁱ_inc are the coherent and incoherent scattering length of ith element and summation index i and j runs over different elements in the unit cell. For Na, σ^Na_coh(= 4πb²_coh) = 1.66b (barns) and σ^Na_inc(= 4πb²_inc) = 1.62b, so both contributions need to be considered. The coherent component is sensitive to two-body correlations, such as in a correlated diffusion process, while the incoherent term reflects the time-correlations of the position of a single atom, or self-diffusion process. Hence, one can quantify the coherent (correlated) and incoherent (independent) contribution of Na diffusion to QENS spectra by simulating S^Na_inc(Q,E) and S^Na–Na_coh(Q,E) given by eqn (13) and (14) from MD:where summation indices i and j run over all Na atoms in the simulation cell.

We computed S^Na_inc(Q,E) and S^Na−Na_coh(Q,E) in Na₃PS₄ with 0.6% Na-vacancy from our NVT-MLMD simulations (Fig. S7, ESI†). Our calculated S^P_inc(Q,E) and S^S_inc(Q,E) do not show any significant quasi-elastic contribution, confirming that Na dynamics dominate the QENS (Fig. S13, ESI†). This can be understood as follows. At elevated temperatures (∼600 K), the large vibrational amplitude of anharmonic low-energy modes (involving Na vibrations and PS₄ librations) results in intermittent stochastic Na hopping along (100) channels, which generates the QENS signal. However, the PS₄ units do not undergo stochastic reorientations at 600 K. Thus, they do not contribute to QENS (from low-frequency stochastic dynamics), but rather to the inelastic signal at higher frequency. To further establish this, we have also calculated an angular autocorrelation function, ς(t) of PS₄ tetrahedra.⁴⁰ This quantity remains very close to unity at all times at 600 K (see Fig. S14, ESI†), confirming the absence of stochastic reorientations of PS₄. The calculated S^Na_inc(Q,E) was fit with a single Lorentzian and the Q-dependence of Γ was satisfyingly described by a CE model. The simulated S^Na–Na_coh(Q,E) has a much smaller intensity than S^Na_inc(Q,E) (Fig. S7, ESI†), except near the Bragg peaks, indicating the dominance of incoherent diffusion. In Fig. 4(c), we show the simulated Γ(Q) from S^Na_inc(Q,E), in excellent agreement with our QENS measurements. Further, we also computed the MSD of each species in Na₃PS₄ with 0.6% Na-vacancy at several temperatures from 100 K to 600 K (Fig. S4, ESI†). The estimated D from MSD is ∼10⁻⁶ cm² s⁻¹ at 600 K, in good agreement with our estimates from QENS. Higher Na-vacancy concentrations enhance Na diffusion in Na₃PS₄ (see Table 3 and Fig. S10, ESI†). A recent study on Na_3−xPn_1−xW_xS₄ (Pn = P, Sb) also showed that Na-vacancy significantly enhances the Na-conductivity,⁶³ consistent with our finding.

Table 3 The calculated Na diffusion constant, D (a) with different degrees of host-flexibility and (b) with different Na-vacancy concentration in Na₃PS₄ at 600 K, on a 2 × 2 × 2 supercell from a linear fit to MSD

(a) Host-flexibility vs. D, with 2% Na-vacancy
Degree of host-flexibility	D (×10^−6 cm^2 s^−1)
Unconstrained	2.4
Frozen PS4	0.5
Frozen PS4 rotation	2.1

---

(b) Na-vacancy vs. D
% Na-Vacancy	D (×10^−6 cm^2 s^−1)
0.0	0.0
0.6	1.0
1.8	2.7

---

So far, we discussed the effect of vacancies and the possible role of low-E phonons on Na diffusion based on our lattice dynamics analysis. We now isolate the effect of host framework dynamics on Na diffusion via constrained MLMD simulations (see Materials and methodology section). We compare NVT-MLMD simulations performed with: (i) unconstrained host-lattice, (ii) frozen host-lattice and, (iii) frozen PS₄ rotation (but ‘wiggling’ translations allowed). The resulting MSD for the three cases are shown in Fig. S11 (ESI†). Upon freezing the host dynamics, D decreases by about a factor of five, clearly demonstrating the importance of lattice flexibility (Table 3). We find that only freezing the PS₄ rotations (case (iii)) also results in a suppression of Na MSD but to a lesser extent than the fully frozen lattice (case (ii)).

Further, we also compare the Na partial DOS in all three cases (Fig. 3(f)). We find that constraining the PS₄ dynamics significantly reduces the soft phonon modes at E ∼ 5 meV. In addition, we find that, conversely, increased Na-vacancy concentrations enhance the phonon-DOS at E ∼ 3–5 meV (Fig. S3, ESI†), and thus the slow PS₄ modes coupled with Na diffusion. Thus, all our simulations point to a significant role of low-energy PS₄ dynamics coupled with Na translation in promoting the Na diffusion in Na₃PS₄. Thus, we have shown that some vacant Na sites are needed to enable diffusion in 1D channels in β-Na₃PS₄. Yet, even in the presence of Na vacancies, the host dynamics of PS₄ units enable a significant increase in the Na diffusivity (see Table 3 and Fig. S11, ESI†). Hence, while a minimum Na-vacancies is required to enable the onset of diffusion, low-energy phonons of the host lattice play an important role of their own.

4 Conclusions

In the present study, we used high-resolution INS and QENS measurements combined with extensive simulations of atomic dynamics, to investigate the role of phonons in enabling fast Na diffusion in Na₃PS₄. We revealed the existence of prominent phonon modes with surprisingly low energy in such a light-atom compound (E < 7 meV), and showed that these modes drastically broaden and soften with increasing temperature, corresponding to a breakdown of harmonic phonon quasiparticles. We computationally established how these modes involve coupled motions of Na ions in 1D channels together with PS₄ dynamics, and how the large population of these modes at elevated temperatures yields a large vibrational amplitude that helps Na to diffuse in 1D channels. Importantly, we show that the dynamics of low energy modes at the Brillouin zone boundary match the minimum energy pathways for Na hopping, and thereby link these anharmonic soft phonon modes with Na diffusion. Further, we investigated the coupled effects of both host dynamics and Na vacancies on diffusion and found that both contribute significantly toward increasing Na diffusivity. Our results emphasize the importance of anharmonic host dynamics for fast diffusion, in addition to the better known effect of vacancies. These results thus represent new insights into the role of phonon dynamics in solid electrolytes. Our QENS measurements and large-scale MLMD simulations both show the presence of long-range diffusion at 600 K and provide diffusion constants in remarkable agreement, D ∼ 10⁻⁶ cm² s⁻¹.

These results establish the importance of anharmonic low-energy modes in a soft crystal structure near a structural phase transition as a means of facilitating fast ion conduction, which provides a new type of descriptor in screening for future solid electrolyte candidates. INS and Raman measurements of phonons as a function of temperature, as well as ultrasonic sound-velocity measurements, can detect soft anharmonic phonon modes and can be used to screen for such effects in candidate SEs.^24–27 In addition, first-principles simulations can be used to reveal incipient lattice instabilities and potentially beneficial anharmonic soft modes.

Conflicts of interest

The authors have no conflicts of interest to declare.

Acknowledgements

We thank Linda Nazar for interesting discussions. Neutron scattering measurements were supported by the USA Department of Energy, Office of Science, Basic Energy Sciences, Materials Sciences and Engineering Division, under Award No. DE-SC0019978. Simulations and data analysis were supported by Duke startup funds. HW acknowledges financial support from NSF EPSCOR RII Track 4 award (No. 2033397). The research at Oak Ridge National Laboratory's Spallation Neutron Source and High Flux Isotope Reactor was sponsored by the Scientific User Facilities Division, Office of Basic Energy Sciences, USA DOE. Theoretical calculations were performed using resources of the National Energy Research Scientific Computing Center, a U.S. DOE Office of Science User Facility supported by the Office of Science of the USA Department of Energy under Contract No. DE-AC02-05CH11231. Z. D. H. was supported by Laboratory Directed Research and Development (LDRD) funding from Argonne National Laboratory, provided by the Director, Office of Science, of the U.S. Department of Energy under Contract No. DE-AC02-06CH11357.

Notes and references

1. J. Janek and W. G. Zeier, Nat. Energy, 2016, 1, 1–4.
2. Z. Zhang, Y. Shao, B. Lotsch, Y.-S. Hu, H. Li, J. Janek, L. F. Nazar, C.-W. Nan, J. Maier and M. Armand, et al. , Energy Environ. Sci., 2018, 11, 1945–1976.
3. X.-B. Cheng, R. Zhang, C.-Z. Zhao and Q. Zhang, Chem. Rev., 2017, 117, 10403–10473.
4. K. H. Park, Q. Bai, D. H. Kim, D. Y. Oh, Y. Zhu, Y. Mo and Y. S. Jung, Adv. Energy Mater., 2018, 8, 1800035.
5. T. Famprikis, P. Canepa, J. A. Dawson, M. S. Islam and C. Masquelier, Nat. Mater., 2019, 18, 1278–1291.
6. A. Manthiram, X. Yu and S. Wang, Nat. Rev. Mater., 2017, 2, 1–16.
7. J. B. Boyce and B. A. Huberman, Phys. Rep., 1979, 51, 189–265.
8. M. D. Slater, D. Kim, E. Lee and C. S. Johnson, Adv. Funct. Mater., 2013, 23, 947–958.
9. N. Tanibata, M. Deguchi, A. Hayashi and M. Tatsumisago, Chem. Mater., 2017, 29, 5232–5238.
10. L. Zhang, K. Yang, J. Mi, L. Lu, L. Zhao, L. Wang, Y. Li and H. Zeng, Adv. Energy Mater., 2015, 5, 1501294.
11. A. Hayashi, K. Noi, N. Tanibata, M. Nagao and M. Tatsumisago, J. Power Sources, 2014, 258, 420–423.
12. 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., 2016, 128, 9786–9790.
13. S. Xiong, Z. Liu, H. Rong, H. Wang, M. McDaniel and H. Chen, Sci. Rep., 2018, 8, 1–7.
14. Z. Zhu, I.-H. Chu, Z. Deng and S. P. Ong, Chem. Mater., 2015, 27, 8318–8325.
15. N. J. de Klerk and M. Wagemaker, Chem. Mater., 2016, 28, 3122–3130.
16. X. He, Y. Zhu and Y. Mo, Nat. Commun., 2017, 8, 1–7.
17. S.-H. Bo, Y. Wang, J. C. Kim, W. D. Richards and G. Ceder, Chem. Mater., 2016, 28, 252–258.
18. M. Duchardt, U. Ruschewitz, S. Adams, S. Dehnen and B. Roling, Angew. Chem., Int. Ed., 2018, 57, 1351–1355.
19. T. Famprikis, J. A. Dawson, F. Fauth, O. Clemens, E. Suard, B. Fleutot, M. Courty, J.-N. Chotard, M. S. Islam and C. Masquelier, ACS Mater. Lett., 2019, 1, 641–646.
20. T. Krauskopf, S. P. Culver and W. G. Zeier, Inorg. Chem., 2018, 57, 4739–4744.
21. C. K. Moon, H.-J. Lee, K. H. Park, H. Kwak, J. W. Heo, K. Choi, H. Yang, M.-S. Kim, S.-T. Hong and J. H. Lee, et al. , ACS Energy Lett., 2018, 3, 2504–2512.
22. Q. Zhang, C. Zhang, Z. D. Hood, M. Chi, C. Liang, N. H. Jalarvo, M. Yu and H. Wang, Chem. Mater., 2020, 32, 2264–2271.
23. T. Famprikis, H. Bouyanfif, P. Canepa, J. Dawson, M. Zbiri, E. Suard, F. Fauth, H. Y. Playford, D. Dambournet, O. Borkiewicz, M. Courty, J.-N. Chotard, S. Islam and C. Masquelier, Chem. Mater., 2021, 33(14), 5652–5667.
24. T. Krauskopf, S. Muy, S. P. Culver, S. Ohno, O. Delaire, Y. Shao-Horn and W. G. Zeier, J. Am. Chem. Soc., 2018, 140, 14464–14473.
25. T. Krauskopf, C. Pompe, M. A. Kraft and W. G. Zeier, Chem. Mater., 2017, 29, 8859–8869.
26. M. A. Kraft, S. P. Culver, M. Calderon, F. Böcher, T. Krauskopf, A. Senyshyn, C. Dietrich, A. Zevalkink, J. Janek and W. G. Zeier, J. Am. Chem. Soc., 2017, 139, 10909–10918.
27. S. Muy, J. C. Bachman, L. Giordano, H.-H. Chang, D. L. Abernathy, D. Bansal, O. Delaire, S. Hori, R. Kanno and F. Maglia, et al. , Energy Environ. Sci., 2018, 11, 850–859.
28. S. Muy, R. Schlem, Y. Shao-Horn and W. G. Zeier, Adv. Energy Mater., 2020, 2002787.
29. S. Ohno, A. Banik, G. F. Dewald, M. A. Kraft, T. Krauskopf, N. Minafra, P. Till, M. Weiss and W. G. Zeier, Prog. Energy, 2020, 2, 022001.
30. K. Gordiz, S. Muy, W. G. Zeier, Y. Shao-Horn and A. Henry, Cell Rep. Phys. Sci., 2021, 2, 100431.
31. K. Wakamura, Solid State Ionics, 2004, 171, 229–235.
32. J. L. Niedziela, D. Bansal, A. F. May, J. Ding, T. Lanigan-Atkins, G. Ehlers, D. L. Abernathy, A. Said and O. Delaire, Nat. Phys., 2019, 15, 73–78.
33. J. Ding, J. L. Niedziela, D. Bansal, J. Wang, X. He, A. F. May, G. Ehlers, D. L. Abernathy, A. Said and A. Alatas, et al. , Proc. Natl. Acad. Sci. U. S. A., 2020, 117, 3930–3937.
34. A. Maradudin, Theory of lattice dynamics in the harmonic approximation, Academic Press, New York, 1971.
35. P. Brüesch, Phonons: Theory and Experiments I - Lattice Dynamics and Models of Interatomic Forces, Springer Series in Solid-State Sciences, New York, 1982. vol. 34.
36. X. He, D. Bansal, B. Winn, S. Chi, L. Boatner and O. Delaire, Phys. Rev. Lett., 2020, 124, 145901.
37. T. Lanigan-Atkins, X. He, M. Krogstad, D. Pajerowski, D. Abernathy, G. N. Xu, Z. Xu, D.-Y. Chung, M. Kanatzidis and S. Rosenkranz, et al. , Nat. Mater., 2021, 1–7.
38. Z. Zhang, E. Ramos, F. Lalère, A. Assoud, K. Kaup, P. Hartman and L. F. Nazar, Energy Environ. Sci., 2018, 11, 87–93.
39. Z. Zhang, P.-N. Roy, H. Li, M. Avdeev and L. F. Nazar, J. Am. Chem. Soc., 2019, 141, 19360–19372.
40. K. Sau and T. Ikeshoji, J. Phys. Chem. C, 2020, 124, 20671–20681.
41. E. P. Ramos, Z. Zhang, A. Assoud, K. Kaup, F. Lalere and L. F. Nazar, Chem. Mater., 2018, 30, 7413–7417.
42. I. Hanghofer, M. Brinek, S. Eisbacher, B. Bitschnau, M. Volck, V. Hennige, I. Hanzu, D. Rettenwander and H. Wilkening, Phys. Chem. Chem. Phys., 2019, 21, 8489–8507.
43. K. E. Kweon, J. B. Varley, P. Shea, N. Adelstein, P. Mehta, T. W. Heo, T. J. Udovic, V. Stavila and B. C. Wood, Chem. Mater., 2017, 29, 9142–9153.
44. J. G. Smith and D. J. Siegel, Nat. Commun., 2020, 11, 1–11.
45. D. L. Abernathy, M. B. Stone, M. Loguillo, M. Lucas, O. Delaire, X. Tang, J. Lin and B. Fultz, Rev. Sci. Instrum., 2012, 83, 015114.
46. O. Arnold, J.-C. Bilheux, J. Borreguero, A. Buts, S. I. Campbell, L. Chapon, M. Doucet, N. Draper, R. F. Leal and M. Gigg, et al. , Nucl. Instrum. Methods Phys. Res., Sect. A, 2014, 764, 156–166.
47. O. Delaire and C. Stassis, Characterization of Materials, ch. Phonon Studies, Wiley, 2012, pp. 1–22.
48. V. F. Sears, Neutron News, 1992, 3, 26–37.
49. E. Mamontov and K. W. Herwig, Rev. Sci. Instrum., 2011, 82, 085109.
50. R. T. Azuah, L. R. Kneller, Y. Qiu, P. L. Tregenna-Piggott, C. M. Brown, J. R. Copley and R. M. Dimeo, J. Res. Natl. Inst. Stand. Technol., 2009, 114, 341.
51. G. Kresse and J. Furthmüller, Comput. Mater. Sci., 1996, 6, 15–50.
52. G. Kresse and D. Joubert, Phys. Rev. B: Condens. Matter Mater. Phys., 1999, 59, 1758.
53. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865.
54. A. Togo and I. Tanaka, Scr. Mater., 2015, 108, 1–5.
55. H. J. Monkhorst and J. D. Pack, Phys. Rev. B: Condens. Matter Mater. Phys., 1976, 13, 5188.
56. T. Tadano, Y. Gohda and S. Tsuneyuki, J. Phys.: Condens. Matter, 2014, 26, 225402.
57. H. Wang, L. Zhang, J. Han and E. Weinan, Comput. Phys. Commun., 2018, 228, 178–184.
58. S. Plimpton, J. Comput. Phys., 1995, 117, 1–19.
59. D. Orobengoa, C. Capillas, M. I. Aroyo and J. M. Perez-Mato, J. Appl. Crystallogr., 2009, 42, 820–833.
60. T. F. Willems, C. H. Rycroft, M. Kazi, J. C. Meza and M. Haranczyk, Microporous Mesoporous Mater., 2012, 149(1), 134–141.
61. X. He, Y. Zhu, A. Epstein and Y. Mo, npj Comput. Mater., 2018, 4, 1–9.
62. M. Bee, Biology and Materials Science, Adam Hilger, Bristol, 1988, p. 193.
63. T. Fuchs, S. P. Culver, P. Till and W. G. Zeier, ACS Energy Lett., 2019, 5, 146–151.

---

Footnotes

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

‡ Present address: Solid State Physics Division, Bhabha Atomic Research Centre, Mumbai 400085, India.

---