Crystal structures, local atomic environments and ion diffusion mechanisms of scandium-substituted NASICON solid electrolytes

The importance of exploring new solid electrolytes for all-solid-state batteries has led to significant interest in NASICON-type materials. Here, the Sc3+-substituted NASICON compositions Na3ScxZr2–x(SiO4)2–x(PO4)1+x (termed N3) and Na2ScyZr2–y(SiO4)1–y(PO4)2+y (termed N2) (x, y = 0–1) are studied as model Na+-ion conducting electrolytes for solid-state batteries. The influence of Sc3+ substitution on the crystal structures and local atomic environments has been characterized by powder X-ray diffraction (XRD) and neutron powder diffraction (NPD), as well as solid-state 23Na, 31P, and 29Si nuclear magnetic resonance (NMR) spectroscopy. A phase transition between 295 and 473 K from monoclinic C2/c to rhombohedral R3c is observed for the N3 compositions, while N2 compositions crystallize in a rhombohedral R3c unit cell in this temperature range. Alternating current (AC) impedance spectroscopy, molecular dynamics (MD), and high temperature 23Na NMR studies are in good agreement, showing that, with a higher ...


INTRODUCTION
Rechargeable lithium ion batteries (LIBs) have come to dominate portable electronics over the past two decades, primarily due to their light weight and high energy density. 1,2 However, there are concerns over the safety of LIBs and the abundance and cost of lithium. Sodium ion batteries (NIBs) are being actively investigated as an alternative, since sodium is more abundant and lower in cost while offering a similar intercalation chemistry to lithium. 3−5 Nevertheless, safety issues related to flammable liquid electrolytes remain a serious concern. In contrast, all solid-state batteries (ASSBs), which use nonflammable ion-conducting solid electrolytes, have been considered as potential candidates for alternative energy storage devices. 6−8 One such family of candidate solid electrolytes is the NASICON (sodium (NA) SuperIonic CONductor) class of materials, which offer both high chemical stability and high Na + mobility within a three-dimensional (3D) framework. 9−17 NASICON-type crystal structures have been widely investigated and debated since the 1980s. 9,10,18−24 The challenges associated with structure determination arise from two main factors: first, the difficulty in obtaining stoichiometric compositions due to possible sodium sublimation at high temperature and corresponding formation of ZrO 2 during synthesis; 25 second, the crystal structure varies according to the synthesis procedure and also to the thermal history of the sample. 26,27 Three models are generally used to describe the 31 P, and 29 Si magic angle spinning (MAS) nuclear magnetic resonance (NMR) spectroscopy investigations together with ab initio magnetic calculations provide insights into both the local structure and Na + ion dynamics. The Na + ion transport properties are then investigated by combined experimental (NMR and alternating current (AC) impedance) and atomistic molecular dynamics (MD) techniques.

METHODS
Synthesis. Na 2 CO 3 , Sc 2 O 3 , ZrOCl 2 ·8H 2 O, SiO 2 (fumed, dried at 773 K for 3 h), and NH 4 H 2 PO 4 were used as raw materials for the synthesis of the target compositions. Na 2 CO 3 , ZrOCl 2 ·8H 2 O, and NH 4 H 2 PO 4 were first mixed in mortar and cold pressed into pellets. The pellets were heated at 443 K for 2 h and then at 673 K for 3 h to evaporate water and decompose NH 4 H 2 PO 4 . After cooling down to room temperature, the pellets were ground, mixed with Sc 2 O 3 and SiO 2 , and finally ball milled (zirconia beads/powder = 3/1 in mass ratio, 400 rpm, 2 h) to yield a homogeneous fine powder. The pristine mixture was cold pressed, fired at 1473 K for 12 h, and pulverized after being cooled down to ambient temperature.
Diffraction. NPD patterns were collected using the high-resolution D1B diffractometer at Institute Laue-Langevin (Grenoble, France). High quality diffraction patterns were recorded between 2θ ranges of 0 to 130°, with a step size of 0.1°, accumulated over 2 h. SXRD experiments at 295 and 473 K were performed at the 11BM beamline (APS, Argonne National Lab, U.S.A.) using a wavelength of 0.414209 Å. Rietveld refinements of the crystal structure including lattice parameters, atomic positions, occupancies and atomic displacement parameters (ADPs), were performed using the Fullprof package. 42 Using the "Bond_Str" interface in the same package, bond valence energy landscape calculations have been performed. In this method, the Rietveld refined structural model is used as input. The unit cell was then considered as a grid containing small units with size of a/10 × b/ 10 × c/10. The energy of locating a Na + ion in each unit was calculated using the bond valence parameters. The bond valence energy landscape map was then obtained by setting an energy threshold (1 eV in our case). Furthermore, the Maximum Entropy Method (MEM)/Rietveld analysis was carried out using the RIETAN-FP software. 43 Differential Scanning Calorimetry (DSC). DSC measurements were performed using a Netzsch DSC 204F1 apparatus operating at heating/cooling rates of 10 K min −1 between 238 and 623 K. During the measurement, samples were heated from 298 to 623 K, cooled down to 238 K, heated to 623 K, and then cooled back to 298 K.
AC Impedance. The ionic conductivity was measured using pellet samples, which had been sintered in a FCT Spark Plasma Sintering apparatus at 100 K min −1 up to 973 K for 30 min under an applied pressure of 100 MPa and under vacuum. The sintered pellets (of 90% compactness) were polished and then metalized on each face by gold sputtering using a Bal-Tec SCD 050. AC impedance measurements were performed under argon atmosphere to avoid any sample contact with oxygen and/or moisture at various stabilized temperatures ranging from 298 to 573 K (upon heating and cooling in steps of 25 K) using a Bio-Logic MTZ-35 Impedance Analyzer. A frequency range of 30 MHz to 0.1 Hz and an excitation voltage of 0.1 V were applied during the measurements. Using a Nyquist plot of the experimental data, the ionic conductivity σ was then calculated using the eq 1: in which h is the thickness of the solid electrolyte pellet, S the pellet surface, and R the total resistance obtained from the measurement. Solid-State NMR Spectroscopy. Powder samples were packed into 2.5, 4.0, or 7.0 mm ZrO 2 rotors (Bruker) and closed with Kel-F or BN caps for ambient temperature (RT) and high temperature (HT) NMR experiments, respectively. HT (320 to 875 K) 23 Na MAS NMR (4 kHz MAS frequency) measurements were performed at 9.4 T (Avance I console) using a Bruker double resonance 7.0 mm MAS Figure 1. Relation between the end members NaZr 2 (PO 4 ) 3 − Na 4 Zr 2 (SiO 4 ) 3 −Na 3 Sc 2 (PO 4 ) 3 and the investigated samples of Na 3 Sc x Zr 2−x (SiO 4 ) 2−x (PO 4 ) 1+x (N3, blue) and Na 2 Sc y Zr 2−y (SiO 4 ) 1−y -(PO 4 ) 2+y (N2, red).

Chemistry of Materials
Article probe with laser heating of the sample. A temperature calibration using KBr was done prior to the measurements. 23 Na NMR signal line shapes were determined by one-pulse experiments with high power pulses of 2.15 μs and a repetition time of 0.1 s. A saturation recovery pulse sequence was applied to determine 23 Na spin−lattice relaxation time constants (T 1 ) at variable temperatures. The 23 Na NMR shifts were referenced to a 1 M NaCl solution in D 2 O. 44 RT 31 P MAS NMR (30 kHz MAS frequency) experiments were performed at 9.4 T using a Bruker 2.5 mm triple resonance MAS probe on an Avance I console. 31 P signals at 0.8 ppm were referenced to ammonium dihydrogen phosphate. 45 A one-pulse sequence with high power 0.75 μs pulses and a recycle delay of 1.0 s was applied to acquire the spectra. RT 29 Si MAS NMR (14 kHz MAS frequency) measurements were performed at 11.7 T (Avance III console) using a Bruker triple resonance 4.0 mm MAS probehead. Signal shifts have been referenced to trimethylsilylpropanoic acid at 0 ppm. 46 The spectra were acquired by applying a one-pulse sequence with the pulse length of 1.30 μs and a recycle delay of 60 s. NMR data acquisition, processing, and raw data handling were done using the Bruker Topspin 3.2 package. Activation energies (E a ) were derived by fits of the temperature dependent 23 Na T 1 values applying both a Bloembergen−Purcell−Pound (BPP, 3D diffusion) and a 2D diffusion model 47−50 using in-house Mathematica scripts. Details on the BPP equation and 2D diffusion model are given in the Supporting Information.
Density Functional Theory (DFT) Calculations of NMR Tensors. Ab initio calculations of local structure and NMR tensors were performed with the CASTEP 8.0 plane wave DFT code. 51 The electron exchange-correlation energy was treated with the Perdew− Burke−Ernzerhof (PBE) 52 generalized gradient approximation (GGA) functional. Valence electrons were described with "on-the-fly" ultrasoft pseudopotentials 53 and Koelling−Harmon scalar relativistic effects were included. Reconstruction of the all-electron wave function with the gauge-including projector augmented-wave (GIPAW) method 54 was carried out to compute NMR tensors in the periodic systems. For all calculations, the Brillouin zone was sampled using a Monkhorst− Pack 55 grid with a k-point spacing of 0.05 × 2π Å −1 and the maximum plane wave kinetic energy cutoff was 700 eV. Convergence was verified with respect to total energy (1 × 10 −3 eV·atom −1 ), 31 P magnetic shielding (±0.1 ppm), and 29 Si magnetic shielding (±0.05 ppm). Chemical shift (δ iso ) and calculated shielding (σ iso ) are related via the linear expression: In principle, if the slope of this line is assumed to be −1, a single comparison between calculation and experiment of a known structure will yield the reference shielding. In practice, when possible, it is more reliable if multiple sites are compared between theory and experiment and linear regression employed to determine the slope and reference shielding. 56 The latter strategy was employed in this work and the relations for both 31 P and 29 Si are taken from the literature; δ iso ( 31 P) = −0.8618 × σ iso ( 31 P) + 237.7 ppm from Diez-Goḿez et al. 57 and δ iso ( 29 Si) = −0.9958 × σ iso ( 29 Si) + 327.3 ppm from a recent work by some of the authors. 58 Starting structures were taken from the refined structural models using synchrotron data. The structures were geometry optimized using DFT forces until the change in energy, maximum atomic force, and maximum atomic displacement were less than 2 × 10 −5 eV·atom −1 , 5 × 10 −2 eV·Å −1 , and 1 × 10 −3 Å, respectively.
Molecular Dynamics Simulations. Ion diffusion at finite temperature was modeled using interatomic potential based molecular dynamics (MD), which has been applied successfully to other battery materials. 49,59−64 The parameters of the interatomic potentials were taken from the extensive library of potentials developed by Pedone et al. 65 and the simulations were performed using the LAMMPS code. 66

Article
Prior to the MD simulations, these parameters have been used to perform energy minimization calculations. The simulated unit cell parameters and bond lengths are close to experimental values with a deviation less than 2%. For each composition, based on the crystal structure determined from the diffraction techniques, a supercell containing about 8,000 atoms (1296 and 864 Na + ions in N3 and N2 type compositions, respectively) was first energy minimized and then MD was applied for 15 ns with a time step of 0.002 ps. Several simulation runs at various temperatures (373 K, 473 K, 573 K, 673 K, 773 K) were performed using the NPT ensemble. To monitor the diffusion, after the system had reached equilibrium, the mean squared displacement (MSD) of Na + ions was resolved as with N equals the number of Na + ions in the simulation cell, r i (t) and r i (0) are the positions of Na + ion number i at time t and time 0, respectively. The chemical diffusion coefficients (D c ) were then calculated from the relation The tracer diffusion coefficients (D t ) were derived in a similar way using the mean squared displacement (MSD). The MSD plots for the sample Na 3 (5) in which f is the correlation factor defined as the ratio of D c over D t , n is the density of charge carrier (number of Na + ions per volume unit), q the charge, k the Boltzmann constant, and T the system temperature. Calculations of Van Hove correlations 67−68 were calculated and plotted according to the procedure given in Supporting Information.

RESULTS AND DISCUSSION
3.1. Crystal Structures and Local Atomic Environm e n t s . D i ff r a ct io n . T h e c r y s t a l s t r u c t u r e s o f Na 3 Sc x Zr 2−x (SiO 4 ) 2−x (PO 4 ) 1+x and Na 2 Sc y Zr 2−y (SiO 4 ) 1−y -(PO 4 ) 2+y (x, y = 0, 0.25, 0.50, 0.75, and 1) have been investigated first using diffraction techniques. For each composition, the synchrotron X-ray diffraction data have been collected at 295 and 473 K. For N3 compositions (Na 3 Sc x Zr 2−x (SiO 4 ) 2−x (PO 4 ) 1+x ) at 295 K, the diffraction patterns can be indexed using the β-NASICON model (C2/c, Supporting Information Table S1). To illustrate this, Figure 2, presents a typical synchrotron XRD pattern of Na 3 Sc 0.25 Zr 1.75 (SiO 4 ) 1.75 (PO 4 ) 1.25 and the Rietveld refinement results. Notably, the pattern cannot be indexed in the higher symmetry γ-NASICON structure (R3̅ c) at 295 K. The reflection at around 2θ = 3.69°is asymmetric and a much more reliable fit is obtained if more than one Bragg reflection is used. The reflection at about 2θ = 7.36 o is not present in the γ model but can be fitted using the β model.
At 473 K, N3 compositions can be indexed in the rhombohedral γ-NASICON model. A reversible thermal phenomenon can be detected in our DSC analysis at about 435 K (Figure 2), confirming the βto γ-phase transition. For N2 compositions (Na 2 Sc y Zr 2−y (SiO 4 ) 1−y (PO 4 ) 2+y ) both the low temperature (295 K) and the high temperature (473 K) diffraction patterns can be indexed using the γ-NASICON model.  3 69 and Na 3 Zr 2 (SiO 4 ) 2 (PO 4 ). 9 One interesting feature of the "NASICON type" compositions is their anisotropic thermal expansion of the unit cell. 19 This can also be observed from the cell parameters listed in Figure 3. When the temperature is increased, minor contraction along the a-and b-directions of the hexagonal cell is observed while the c-axis undergoes significant expansion. The net result is that the value of V/Z does not change much with temperature. In the most significant case, a 1.2% increase is observed from 295 to 473 K. Such a small volume expansion is an important criterion for all-solid-state batteries, since large changes in volume across the solid−solid interfaces can lead to delamination which breaks the sodium ion pathways.
In the NASICON crystal structure, Sc 3+ and Zr 4+ share the same crystallographic site and form MO 6 octahedra, and Si 4+ and P 5+ share the same crystallographic site and form XO 4 tetrahedra. Two MO 6 octahedra are corner-connected to three XO 4 tetrahedra, forming the well-known "lantern units" ( Figure  4). In the β-NASICON crystal structure, three crystallographically different Na sites are present. In contrast, in the γ-NASICON phase, because of the symmetry there are only two types of Na sites, Na(1) and Na (2). In our refinement of the high temperature (473 K) structure, when the ADPs are refined anisotropically, the Na(1) site forms a flat disc in the ab-plane while Na(2) site forms a long ellipsoid along the c-direction. Their equivalent isotropic ADP values are listed in Figure 5 and in Supporting Information Table S3. As x or y increases, a   23 Na NMR. Room temperature (RT) 23 Na MAS NMR spectra show one symmetric signal at approximately −28 ppm with a number of rotational sidebands (RSBs) ( Figure 6). The Na local environments are too chemically similar to be resolved at the applied MAS speed of 4 kHz and under the moderate magnetic field strength of 9.4 T. The residual line broadening coming at least in part from the second-order quadrupolar interaction is not completely removed by MAS. The 23 Na NMR signal line shape, shift, and full width at half-maximum (fwhm) change significantly with increasing temperature (Figure 7). When the temperature is increased above 374 K, the RSBs are no longer observed. The fwhm reaches a maximum at 374 K, dropping steadily thereafter until it reaches a plateau above 500 K. The average chemical shift moves steadily to lower frequencies until a temperature of 500 K is reached; thereafter the shift remains essentially constant.

Chemistry of Materials
The high temperature (HT) 23 Na NMR line shape development indicates a strong influence of the temperature on the Na + ion diffusion. The increase of fwhm between RT Key: ScO 6 /ZrO 6 , gray octahedra; SiO 4 /PO 4 , red tetrahedra; Na1, yellow discs; Na2, blue ellipsoids. Na1 and Na2 are shown as displacement ellipsoids with a probability value of 95%.

Chemistry of Materials
Article and 374 K is ascribed to either a change in the Na local environments of sites occupied by Na with increasing temperature or to an interference between motion and magic angle spinning motion on the MAS time scale. This prevents the method from removing sources of line broadening such as the quadrupolar interaction and dipolar coupling. 70 The steady decrease in fwhm above 374 K indicates increased motion. The narrowing limit (or fast regime) is reached above 500 K where ionic motion is occurring on a time scale that effectively averages any shift distribution, dipolar coupling, and quadrupolar coupling interactions, i.e., the hops between different Na + sites are expected to occur with frequencies in the kHz regime or above, on the basis of typical values for these interactions. 31 P NMR. Room temperature 31 P MAS NMR reveals two signals for Na 2 Zr 2 (SiO 4 )(PO 4 ) 2 and Na 3 Zr 2 (SiO 4 ) 2 (PO 4 ) (Figure 8) of which the minor intensity resonance peaks are due to the formation of a zirconium-deficient impurity phase that has been reported previously. 71 Only one asymmetric 31 P signal is present for all other compositions ( Figure 8). One might expect the asymmetry of these signals to be caused by the random distribution of Zr/Sc on the second coordination spheres with respect to the P site, as recently shown in the NASICON structure Li(Ti/ Ge) 2 (PO 4 ) 3 solid solution. 57 This Zr/Sc disorder leads to various local chemical environments: P(OZr) 4 , P(OZr) 3 (OSc), P(OZr) 2 (OSc) 2 , P(OZr)(OSc) 3 , and P(OSc) 4 . Ab initio calculations were undertaken to understand the local structure and NMR results, which revealed that the 31 P chemical shift is more strongly affected by the next-nearest neighbor (NNN) sodium occupancy than the NNN metal distribution (Supporting Information Figure S4). Several spectral features may be understood in this context. First, the fwhm of the 31 P NMR resonances of Zr end-members (x, y = 0.00) Na 2 Zr 2 (SiO 4 )(PO 4 ) 2 and Na 3 Zr 2 (SiO 4 ) 2 (PO 4 ) were 5.5 and 3.5 ppm, respectively, significantly larger than the fwhm of ≤ 0.5 ppm observed for LiTi 2 (PO 4 ) 3 and LiGe 2 (PO 4 ) 3 in the fully lithium-ordered NASICON phases in prior studies. 57,72 Furthermore, these 31 P resonances are asymmetric and contain shoulders, which cannot be explained by Zr/Sc disorder since Sc is not present at x, y = 0 but can be rationalized by the Na disorder. The ab initio NMR calculations from relaxed structures yielded a range of 31 P shifts from −11 to −25 ppm in Na 2 Sc y Zr 2−y (SiO 4 ) 1−y (PO 4 ) 2+y and −7 to −19 ppm for Na 3 Sc x Zr 2−x (SiO 4 ) 2−x (PO 4 ) 1+x , in excellent agreement with experimental results. The slight positive shift from N2 to N3 is also explained by the average increase in Na NNN coordination in N3; approximately, for both series, each sodium NNN (defined by a 4.5 Å radius) shifts the 31 P resonance by approximately +4 ppm. 29 Si NMR. Room temperature 29 Si MAS NMR reveals two signals at −94 and −111 ppm (Supporting Information Figure  S5). The resonance at −111 ppm is due to SiO 2 impurity, which could not be detected in the diffraction experiments due to its amorphous nature. The −94 ppm signal intensity continuously decreases with higher x values for Na 3 Sc x Zr 2 − x (SiO 4 ) 2 − x (PO 4 ) 1 + x and y value s for Na 2 Sc y Zr 2−y (SiO 4 ) 1−y (PO 4 ) 2+y , respectively, which is in agreement with the substitution of the SiO 4 /PO 4 units. Ab initio magnetic shielding calculations suggest a range of −93 to −101 ppm in Na 2 Sc y Zr 2−y (SiO 4 ) 1−y (PO 4 ) 2+y and −85 to −97 ppm for Na 3 Sc x Zr 2−x (SiO 4 ) 2−x (PO 4 ) 1+x for the 29 Si NMR spectra, which supports the experimental assignment and observed range. These calculations confirm a silicon shift dependence analogous to that of phosphorus; each sodium NNN results in a shift of the 29 Si resonance by approximately +4 ppm. Thus, unlike the lithium-ordered Li(Ti/Ge) 2 (PO 4 ) 3 solid solution, 57 the cation distribution may be a contributing factor to the chemical shift; however, it is relatively modest in our samples and cannot be distinguished in the sodium disorder-dominated NMR spectra.
3.2. Na + Ion Conductivity. Na + diffusion rates and conductivities are important properties for the effective operation of a solid state battery, which we investigate here using AC impedance, MD simulations, and 23 Na NMR spin− lattice relaxation time measurements. The ionic conductivity has been measured using AC impedance for three different compositions: Na 3 Sc 0.25 Zr 1.75 (SiO 4 ) 1.75 (PO 4 ) 1.25 , Na 2 Zr 2 (SiO 4 )-(PO 4 ) 2 , and Na 2 ScZr(PO 4 ) 3 . Figure 9 shows typical Nyquist plots obtained for Na 3 Sc 0.25 Zr 1.75 (SiO 4 ) 1.75 (PO 4 ) 1.25 at various temperatures. An elongated semicircle is observed at high frequencies of less than 1 kHz, which can be attributed to Naion conduction, and an oblique line at lower frequency due to the ionic blocking (gold) electrode. The spectra were fitted using the equivalent circuit inserted in Figure 9a.
At low temperatures and high frequencies, the contribution from the bulk of the material and the grain boundary can be separated (corresponding to R2//CPE2 and R3//CPE3, respectively). This separation can be observed up to 373 K. A fast reduction of the semicircle diameter is noticed upon heating, revealing the enhancement of the ionic conductivity  with temperature. The measured capacitance gives values of between 10 −10 and 10 −9 F, in agreement with the value suggested by Irvine et al. 73 The evolution of the ionic conductivity at various temperatures is plotted in Figure 9b (sum of bulk and grain boundaries contributions) for Na 3 Sc 0.25 Zr 1.75 (SiO 4 ) 1.75 (PO 4 ) 1.25 . We note the superposition of ionic conductivity values between heating and cooling steps, consistent with a process that is reversible with temperature, and that there are no additional detrimental reactions during the measurements.

Chemistry of Materials
From these conductivity values, activation energies for Na +ion diffusion can be obtained using an Arrhenius plot (log(σT) vs 1/T), which is shown in Figure 10. An inflection point can be observed in the Arrhenius plot around 448 K for  (Figure 10, red), indicating a modification of the activation energy. This can be linked to the phase transition from monoclinic C2/c to rhombohedral R3̅ c, previously observed for this type of compound in which the phase transition induces a modification of ionic transport properties. No inflection point in the Arrhenius plot is observed for N2 compounds (Figure 10, blue and green), which is in agreement with the results from diffraction. For N2 compounds, the addition of Sc 3+ increases the activation energy from 0.35 to 0.50 eV and decreases the ionic conductivity. The influence of Sc 3+ -substitution on the ionic conductivity is further revealed by the MD simulated ionic conductivities at 473 K, summarized in Table 1.
From these ionic conductivity values, two key features appear. First, with the same Sc 3+ content per formula unit (x = y), N3 has a higher ionic conductivity than N2 (except for the case of x = y = 1, where the two conductivity values are very close to each other); this result is probably due to the increased number of ionic charge carriers in the N3 system. Second, the ionic conductivity decreases as more Sc 3+ is substituted on the Zr 4+ site, in agreement with the AC impedance measurements. The maximum conductivity is achieved in the composition Na 3 Zr 2 (SiO 4 ) 2 (PO 4 ) (N3 x = 0, with a Si/P ratio of 2/1). The trend is in line with that of the isotropic ADP values discussed in section 3.1. With a higher x or y value there are two key features: (1) the decrease of the lattice parameters leads to a

Chemistry of Materials
Article more compact unit cell, leaving less free volume for sodium-ion diffusion; (2) the crystal structure is more disordered due to the random distribution of Sc/Zr on M(O 6 ) sites and of Si/P on X(O 4 ) sites. This suggests that the crystal structural disorder impedes sodium-ion diffusivity. Roy et al. 61,74 have shown in their MD simulations that the ordering of XO 4 groups significantly influences the Na + conductivity: depending on the distribution of XO 4 groups, the ionic conductivity can vary by up to a factor of 40.
A poor fit to the 23 Na NMR spin−lattice relaxation times (T 1 ) is obtained with a BPP model (3D diffusion), particularly in the HT regime ( Figure 11). On the other hand, the application of a 2D diffusion model results in a very good agreement of NMR data and fit. This indicates that the nature of the Na + motion giving rise to the T 1 temperature dependence with a time scale of 10 −8 s at approximately 416 K (T 1 minimum 47 ) is 2D on the local atomic scale of NMR although the long-range diffusion pathways are 3D. This finding is in agreement with the rapid in-plane motion of Na + ions at the Na1 sites, as revealed by the disc-shaped ADP (Figure 4). A detailed discussion of the dimensionality of Na + diffusion can be found in the Supporting Information.
The activation energies (E a ) have been derived from NMR, AC impedance, and MD simulations, which are summarized in Supporting Information Table S7 and Figure S8. All E a values are between 0.2 and 0.5 eV, which are typical values for solid electrolyte materials. [6][7][8]36,75 In general, the results from all the techniques demonstrate the same trends across the composition range: higher activation energy values are obtained for larger x (or y) values (lower Si/P ratios). This agrees well with the trends in the ionic conductivity and confirms the good agreement between experiments and simulations. In the recent work of Samiee et al., DFT calculations confirmed that Na +motion activation barriers increase with a decrease in Si/P ratio in the monoclinic NASICON structure. 76 The diffusion coefficients D T obtained from NMR and the chemical diffusion coefficients D c calculated from the MD simulations have similar trends and values as a function of composition. For all compositions, the diffusion coefficient has a value of the order of 10 −8 cm 2 /s at 473 K, as shown in Figure  12. The diffusion coefficient decreases with a higher concentration of Sc 3+ -substitution, showing the same trends with that of the ionic conductivity data.

Na + Ion Diffusion Mechanism. Diffusion Pathways.
The Na + diffusion mechanism has been studied by both experimental and computational techniques to investigate the Na + diffusion pathways. 10,21,61,74 In Figure 13, we plot the conduction pathways of Na + ions in Na 3 Sc 0.25 Zr 1.75 (SiO 4 ) 1.75 -(PO 4 ) 1.25 using three different approaches: (i) MD density plots of the accumulated sodium ion trajectories over the simulated time scale, which reveal the migration pathways and the regions in the lattice that are most frequently traversed by the mobile Na + ions. (ii) Bond valence energy landscapes (BVEL) using the bond valence sum method, 77,78 which has been applied to a range of ionic conductors. 79,80 This method requires only the crystal structure as an input and can probe likely ion diffusion pathways with minimal computational cost, but lattice relaxation effects or cooperative mechanisms are not taken into account. (iii) The maximum-entropy method (MEM)/Rietveld analysis uses the scattering densities by giving the maximum variance of calculated structure factors within standard deviations of observed ones. In this way, structural information can be extracted effectively from the    Figure 13). Such 3D diffusion behavior of sodium ions allows access to all possible pathways and helps to enhance the ionic conductivity. Furthermore, the Na + density plot forms flat discs on the Na(1) sites and thin long ellipsoids on the Na(2) sites, which is exactly the same feature that is present in the Rietveld refined anisotropic atomic displacement parameters.

Chemistry of Materials
Na(1)−Na(2) Site Exchange. As mentioned in Section 3.1, the rhombohedral NASICON form has two crystallographically different Na sites. The local Na−Na environment is shown in Figure 14, where each Na(1) site is surrounded by six Na (2) sites and each Na(2) site is surrounded by two Na(1) sites. The approximate distances between different sites are indicated.
Given that the Na(1) sites have a large separation (6.48 Å) and the Na(2) sites are also far apart (4.98 Å), it is reasonable to consider that most Na migration occurs between Na(1) and Na (2) sites where the distance is 3.46 Å. To verify such site exchange we re-examined the MD density plots. In Figure 13 the atomistic trajectory has been recorded separately for Na +ions according to their initial site: the density plots of Na + ions initially located at Na(1) and Na (2) sites at the beginning of the MD simulation are presented on the left and right of Figure  15, respectively. As shown, when the two density plots are overlaid on top of each other, they are very similar and form well connected 3D pathways between the Na(1) and Na (2) sites. Furthermore, the zigzag pathways are in good agreement with the result reported by Kabbour et al. 29 using probability density function calculations.
Na + -Ion Intersite Hopping. The typical time scales of sodium intersite hopping have not been extensively investigated previously for NASICON materials. To extract this information from the MD results we have calculated the van Hove correlation functions of which the transformed version of the self-part r 2 G s (r, t) is related to the probability of an atom traveling a distance r after a time interval of t. In Figure 16 we show an example of the r 2 G s (r, t) of Na 3  Four main local maxima can be observed, from low to high r values: (i) A first sharp and intense peak appears at r = 0.5 Å which is related to local on-site vibrations of the Na + ions. A slight right shift of this peak is observed as the temperature of the simulation system increases, in line with the increased

Chemistry of Materials
Article vibrational amplitudes as the kinetic energy of the Na + ions is raised. (ii) The second peak appearing at about r = 3.5 Å corresponds to the distance between the Na(1) site and the Na(2) site as shown in Figure 14. This peak is indicative of hopping between neighboring Na(1) and Na(2) sites in accord with the NMR results. (iii) The third peak at about r = 5 Å is the distance between two neighboring Na(2) sites. As there is an oxygen atom located directly between two Na(2) sites, direct migration of Na + between Na(2) sites is unfavorable. We observe instead that the Na + ion goes through a two-step process: starting from an Na(2) site, it moves to an Na(1) and then continues to another Na(2) site. (iv) The next peak at r = 6.5 Å is the distance between two Na(1) sites. Similar to the previous case, the peak can be considered as a two-step migration process following sites Na(1)−Na(2)−Na (1).
By plotting the self-part of the van Hove correlation functions with different time interval values ( Figure 17) we can obtain insights into the time dependence of the Na + ion diffusion process. The van Hove correlation function starts to form a peak at 3.5 Å in the first 10 ps, indicating that the site exchange between Na(1) and Na(2) sites can occur within a very short time scale. This peak has its highest intensity on time scales of t = 250 ps. As the time interval value increases, Na + ions are able to accomplish the two-step migration process (Na(1) → Na(2) → Na(1) or Na(2) → Na(1) → Na(2)), resulting in a decrease of intensity of the peak at 3.5 Å but an increase of intensity of the peaks at 5 and 6.5 Å. The overall picture is of short-range hopping between Na(1) and Na(2) sites on a time scale of a few hundred picoseconds (where some of this is hopping back to the initial site) and long-range hopping on time scales of nanoseconds.

CONCLUSIONS
The NASICON-type Na 3 Sc x Zr 2−x (SiO 4 ) 2−x (PO 4 ) 1+x and Na 2 Sc y Zr 2−y (SiO 4 ) 1−y (PO 4 ) 2+y solid electrolytes have been investigated through a powerful multitechnique approach to obtain detailed information on the crystal structures, local atomic environments, and Na + ion diffusion mechanisms.
First, the bulk crystal and local structures of these compositions have been determined using a combination of diffraction and solid-state NMR techniques. A phase transition between 295 and 473 K from monoclinic C2/c to rhombohedral R3̅ c is observed for the Na 3 Sc x Zr 2−x (SiO 4 ) 2−x -(PO 4 ) 1+x compositions, while Na 2 Sc y Zr 2−y (SiO 4 ) 1−y (PO 4 ) 2+y materials crystallize in a rhombohedral R3̅ c unit cell in the same temperature range. Second, both AC impedance measurements and molecular dynamics (MD) simulations show ionic conductivity values of 10 −5 S/cm at 298 K and 10 −3 S/cm at 473 K for these Sc 3+ -substituted NASICON materials. Finally, the 3D diffusion pathways of sodium ions in the NASICON structure have been revealed, for the first time, by three different approaches: MD, bond valence, and maximum entropy methods. 23 Na NMR measurements reveal that the 3D diffusion pathway consists of several 2D local hops between Na(1) and Na(2) sites. A key feature is that ion diffusion proceeds by site exchange between Na(1) and Na(2) sites with long-range diffusion occurring on time scales of less than a nanosecond (at 773 K).
This systematic study provides important insights into the wider family of NASICON-type materials and will be valuable in identifying next-generation solid electrolytes.

* S Supporting Information
The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.chemmater.7b05237.
Unit cell parameters, GoF, Rwp, Rp, and ADPs values of the Rietveld refinements; 29 Si NMR signals from measurements and DFT calculations; MSD plots, diffusion coefficient, ionic conductivity, and activation energy values from MD simulations; models used for NMR; and discussion of the dimensionality of the Na + diffusion (PDF) Figure 16. Na + -ion dynamics in Na 3 Sc 0.25 Zr 1.75 (SiO 4 ) 1.75 (PO 4 ) 1.25 simulated at 773 K using a time scale of 500 ps. Vertical axis is the transformed self-part of the van Hove correlation function r 2 G s (r, t) which is a measure of probability of moving a distance r. Different Na + dynamical processes are indicated, and the colors of these correspond to the site exchange process indicated in Figure 14.

Chemistry of Materials
Article ■ AUTHOR INFORMATION