Theoretical and Experimental Studies of ion Transport in Mixed Polyanion Solid Electrolytes

and display some of the highest ionic conductivities reported to date. However, the effect of polyanion mixing on ion transport properties is still debated. Here, we focus on Na 1+x Zr 2 Si x P 3−x O 12 (0 ≤ x ≤ 3) NASICON electrolyte to elucidate the role of polyanion mixing on Na-transport properties. Although there is a large body of data available on this NASICON system, transport properties extracted from experiments or theory vary by orders of magnitude, signifying the need to bridge the gap between different studies. Here, more than 2,000 distinct ab initio -based kinetic Monte Carlo simulations have been used to map the statistically vast compositional space of NASICON over an unprecedented time range and spatial resolution and across a range of temperatures. We performed impedance spectroscopy of samples with varying Na compositions revealing that the highest ionic conductivity (~ 0.1 S cm –1 ) is achieved in Na 3.4 Zr 2 Si 2.4 P 0.6 O 12 , in line with our predictions (~ 0.2 S cm –1 ). Our predictions indicate that suitably doped NASICON compositions, especially with high silicon content, can achieve high Na + mobilities. Our findings are relevant for the optimization of mixed polyanion solid electrolytes and electrodes, including sulfide-based polyanion frameworks, which are known for their superior ionic conductivities.


Introduction
The reliability of rechargeable batteries and other energy storage devices depends on the fast delivery of ions between electrodes. [1][2][3][4] In rechargeable batteries, the safety and performance of electrolytes are as important as other battery components, such as electrodes. Solid-state batteries, where liquid electrolytes are replaced by solid fastion conductors, offer a promising pathway for safer commercial lithium-ion and prototypical sodium (Na)-ion batteries. [4][5][6] Several solid electrolytes, such as the LiSiCON-type Li4±xSi1-xZxO4 (Z = P 5+ , Al 3+ , Sn 4+ and/or Ge 4+ ), 7,8 Li10MP2S12 made from Li4MS4:Li3PS4 (M = Ge 4+ , Sn 4+ , and Si 4+ ), [9][10][11][12][13][14][15][16] and Na1+xZr2SixP3−xO12 (hereafter referred as NASICON), [17][18][19] leverage mixed polyanion frameworks to attain some of the highest ionic conductivities (≈10 mS/cm) reported so far. 4 Ion transport in mixed polyanion solid electrolytes has been investigated using both computation and experiments before. 4,9,18,7,11 Yet transport properties measured from different techniques typically vary by orders of magnitude. 20 On one hand, averaged transport properties obtained from experiments, such as impedance spectroscopy, solid-state nuclear magnetic resonance, and quasi-elastic neutron scattering, 4,6,16,21 may incorporate effects arising from defects and grain boundaries. 22 On the other hand, classical or ab initio molecular dynamics (MD) studies perform "one shot" simulations on selected bulk structures, which may not be sufficiently representative in terms of longer time-scale transport processes as well as in the number of possible arrangements of different polyanions. 2,3,[23][24][25] Thus, there is a need to reconcile the measurements and simulations to guide the development and manufacturing of the next generation of solid electrolytes and their batteries.
In this study, we have selected the NASICON Na1+xZr2SixP3−xO12 (0 ≤ x ≤ 3) electrolyte as an example to elucidate the role of polyanion mixing on the macroscopic transport properties, including ionic diffusivity and conductivity. The choice of NASICON is justified by the large body of data available on this system since its discovery more than 40 years ago, 17,18,[26][27][28][29] making it easier to reconcile previous experimental and computational studies.
To capture the large statistical variance in transport properties introduced by mixed polyanions in NASICON samples, we developed a high-fidelity kinetic Monte Carlo (kMC) model bearing the accuracy of density functional theory (DFT) calculations. More than 2,000 distinct kMC simulations served to map the statistically vast compositional space of Na1+xZr2SixP3−xO12 over an unprecedented time rangein the realm of milliseconds -and spatial resolution, with varying temperature.
First of all, our model reproduces existing measurements of Na-transport properties in NASICON, suggesting that a robust sampling of both the spatial and temporal axes is required in computations to accurately estimate these properties. The reproduction of measured Na-transport properties in NASICON also implies that our model correctly captures the collective nature of Na-transport, which is responsible for the high ionic conductivity observed in NASICON. Second, our simulations elucidate the impact of the thermodynamic forces driving the random distribution of PO4 3and SiO4 4groups in NASICON during synthesis (and subsequent thermal treatment) on the Na-ion transport. Third, our statistical insights can guide the selection of optimal doping and thermal treatment strategies to further improve the properties of NASICON electrolytes. For example, our analysis suggests that higher Na + conductivity in NASICON can be achieved by increasing the content of SiO4 4units in place of PO4 3-moieties, while maintaining high Na content in the structure. Motivated by our computational findings, we have also synthesized and characterized selected compositions of the NASICON, which validates our predictions. These findings, elaborated on NASICON electrolytes, are general and transferable to the study of ionic transport in other topical mixed-polyanion Li and Na solid electrolytes. Figure 1a depicts the crystal structure of NASICON, where red SiO4/PO4 tetrahedra share corners with two blue ZrO6 octahedra forming the so-called "lantern units". In the rhombohedral high-temperature NASICON structure ( 3 � ), there are two independent and partially occupied sodium sites, Na(1) and Na (2). The partial occupancy of these sites gives rise to the well-known high ionic conductivity in NASICON. A third interstitial Na site, known as Na (3), adjacent to Na(1) and Na (2) has been previously considered. 30,31 However, Zhang et al. 32 showed that Na (3) (1) sites are indicated by silver spheres, the Na(2) by orange spheres, the (Si/P)O 4 groups by red tetrahedra, Si/P atoms by red spheres, and ZrO 6 units by blue octahedra. b and c depict the local environment of Na(1), with each Na(1) surrounded by six neighboring Na(2) atoms (orange spheres) and six Si/P (red spheres) atoms. For simplicity, O and Zr atoms are not shown in b, c, and d. Each silver hexagonal prism in b or c represents the first coordination shell of a Na(1) site. Panel d is the migration unit used to study Na-migration in NASICON, and Na must hop across several different migration units to ensure Na diffusion. Red triangles in d indicate the bottlenecks caused by Si/PO 4 tetrahedra (oxygen atoms are not shown). e shows the averaged kinetically resolved activation (KRA) barriers for Na(2)⟷Na(1) hops, with varying Na(2) site occupation and Si/P content per migration unit. The barriers were extracted from a local cluster expansion model, which was fitted to the calculated nudged elastic band (NEB) barriers. The diagonal line in e shows locally charge neutral Si/P configurations. The computed NEB barriers used for fitting to generate the heatmap in e are available in the supporting information (SI).

Figure 1d
shows the smallest relevant unit -the migration unit-encompassing all possible migration events for Na-ions in NASICON. The migration unit centered around Na(1) highlights the connectivity between the Na(1) site and its closest six Na(2) neighbor sites, as well as the neighboring Si(P) atoms. Each migration unit (see Figure 1b and c) is edge-shared with six neighboring migration units, where each unit can exhibit its own distribution of Na/vacancies and Si/P. In the NASICON structure, Na-transport must always involve both Na(1) and Na(2) sites, 18,30,33,34 giving rise to migration pathways of the kind Na(1) ⟷ Na (2). The specific Na-migration pathway depends entirely on the Na and Si/P content within the structure, 19,30,31 which sets the occupation and the relative stability of Na(1) and Na(2) sites. For example, at low Na content (x = 0), Na ions reside only in Na(1) sites, favoring a Na(1) ⟷ Na(2) ⟷ Na (1) pathway. 17,30,32 At higher Na content (x ≥ 2), both Na(1) and Na(2) sites are occupied to varying degrees (see later Figure 2), and any of the two migration pathways can be active. 19,30 Importantly, for macroscopic Na diffusion to occur in NASICON, several successful migration hopping events, between adjacent migration units, need to happen, for which it is critical that Na atoms can hop from Na(1) to Na(2) sites (or viceversa) within each unit. Thus, we considered a single migration event to be Na(1) ⟷ Na(2) within a single migration unit, in the presence of varying Na and Si/P content.
For any Na composition in the range 0 ≤ x ≤ 3, the concentrations of Si, and P are set by the NASICON stoichiometry, such that the composition guarantees global charge neutrality (i.e., the NASICON cell is electrostatically neutral). However, each migration unit can exhibit random occupancies of Si/P in Si/P sites, and Na/vacancies in the Na(2) sites (panels b, c and d of Figure 1). Consequently, the migration of Naions may occur through migration units that are not locally charge neutral, i.e., a given unit by itself may not be electrostatically neutral. In our work, we consider all possible Na(1)↔Na(2) hopping events within a given migration unit, including scenarios with and without local charge neutrality.
We evaluated all the unique Na-ion migration pathways that are possible within a migration unit of Figure 1d using the nudged elastic band method (NEB) 35 together with density functional theory (DFT), and using the strongly constrained and appropriately normed (SCAN) meta-generalized gradient approximation. 36 The NEB method is designed to calculate migration barriers in solid electrolytes, as highlighted by prior studies. 2,3,37,38 We identified 32 unique Na-hops, given the migration unit of Figure 1d, which included barriers computed at different Si and P compositions and Na/vacancy contents (Section S2 of the Supporting Information, SI). The calculated migration barriers at 0 K are defined as the difference between the highest energy state and the configuration with the lowest energy --the initial and/or final end point. 37 The accuracy of our computed barriers is validated by the excellent agreement with existing data. For example, our calculated barrier at NaZr2(PO4)3 (~ 458 meV, Table   S1) agrees well with existing experimental values (~ 470 meV) 26 and ab initio MD (~ 472 meV). 32 Typically, migration events have an intrinsic directional dependence if the initial and final states are not equivalent, i.e., a Na(1) → Na(2) hop may exhibit a different migration barrier compared to the reverse hop, Na(2) → Na(1), due to energy differences between the end points, namely, Na(1) and Na (2). Van der Ven et al. 39 introduced the kinetically resolved activation (KRA) barrier, or E (see Equation 2 in the Methods section) that removes the directional dependence of the migration (and migration barriers) with Si content can be attributed to lower Si 4+ -Na + electrostatic repulsion compared to P 5+ -Na + repulsion during Na-migration.
Changing the Na concentration at a given Si/P content within the migration unit does not result in a monotonic increase or decrease in average E . For instance, at 4/6 Si @ Si/P site in Figure 1e, the average E is low at both 0/6 Na @ Na(2) site and at 6/6 Na @ Na(2) site, with the average E reaching a maximum at 3/6 Na @ Na (2) site. This non-monotonic behavior is likely due to a combination of local electrostatic repulsions and local charge imbalance within the migration unit. For example, at 6/6 Na @ Na(2) site and 0/6 Si @ Si/P site (bottom right corner of Figure   1e), a high E can be attributed to the higher electrostatic repulsion between the Na + and P 5+ ions. In contrast, low barriers at both 0/6 Na @ Na(2) and 6/6 Na @ Na (2) along 4/6 Si @ Si/P are likely because of high local charge imbalance that destabilizes the initial/final Na(1) and Na (2) configurations along the pathway. Eventually, high Na content in the migration unit contributes to low E , particularly at high Si content, as indicated by the blue squares towards the top right corner of Figure 1e. Indeed, the lowest values of E are observed at both high Na and high Si contents.
Combining the LCE and kMC, we investigated the relevant Na-transport properties in NASICON, such as ionic diffusivity and conductivity. The main advantage of using a kMC framework is the access to statistically significant amounts of data over large length and time scales, as well as the ability to sample a wide variety of configurations. By comparison, a classical or ab initio MD simulation 3,23-25 performs a "one shot" representative calculation of a system at a given atomic configuration and composition, over smaller length and orders of magnitude shorter time scales.
In total, more than 2,000 kMC simulations were performed to predict the Na + diffusivities and conductivities at selected temperatures of 373, 473 and 573 K of 50 NASICON models (with varying Si/P configurations) at 11 distinct Na concentrations.
Using canonical Monte Carlo simulations from our previous work, 19 550 starting NASICON configurations (50 configurations × 11 concentrations) were generated at high temperature (~1500 K), which mimics the typical synthesis conditions and thermal where is the elemental charge, is the concentration of Na-ions, B is the Boltzmann constant and is the temperature. Importantly, in Equation 1, D is the diffusivity of the center-of-mass of Na-ions in the NASICON structure, which includes cross-correlations of the migrations of different Na-ions (see Section S4 of SI) and not the tracer diffusivity, as often used in many studies. 23,39,40 As confirmed by Figure 2a and Figure 2b, a temperature increase can significantly boost both Na + diffusivities and conductivities across the full Na concentration. Interestingly, Na diffusivity and conductivity increase for increasing Na (and Si) concentrations up to x ~ 2.4, beyond which there is a drop in both diffusivity and conductivity. Thus, we conclude that Na-ion transport in Na1+xZr2SixP3−xO12 is fastest at x ~ 2.4, at all temperatures considered in this work, which agrees with our experimental observations presented below. Figure S6 reveals a significant spread of experimentally measured ionic conductivities from different reports, 18,27,33,[42][43][44][45] which is indicative of the difficulties in performing accurate measurements as well as the large variability in sample preparation. Nonetheless, for Na3Zr2Si2P1O12 our measured experimental conductivity (see Figure 4b) of ~ 0.1 S cm -1 at 573 K is in excellent agreement with our computed value of ~ 0.2 S cm -1 . Our model also agrees with previous measurements 18,19,30,32 which report a low conductivity and diffusivity in Na1Zr2P3O12.
Previous studies have shown that correlated ion-migration can contribute to high diffusivity and conductivity in ionic conductors. 3,4,7,23,32,40,46,47 In regimes of nondilute diffusion carriers, the Haven's ratio (HR) of Figure 2c quantifies the degree of cross-correlation between migrating Na + ions, i.e., the extent to which individual Na + hops contribute to the overall motion of the Na + center-of-mass (Equation S8 of SI).
HR varies between 1 (no cross-correlation or a random movement of individual Na + ions) and 0 (fully correlated transport or non-random movement  (1) and Na (2), 30 which are separated by shorter distances (~ 3.48 Å) than Na(1)-Na (1), possibly facilitating correlated migration of Na + ions. Interestingly, cross-correlation between Na + ions shows an increase with increasing temperature, as indicated by the drop in HR with rising temperature at all x (Figure 2c); this trend also coincides with increasing of D and σ, highlighting the positive impact of correlated motion on the overall Na-ion transport.
Another type of correlation, specifically between successive jumps of the same Na + ion, is to estimate the average correlation factor, f, of Figure 2d (see Equation

S10
of the SI). f measures the variation in local environments (resulting in a change of migration barriers) as experienced by individual Na + ions as they migrate through NASICON. f captures the deviation away from a truly random walk of a given Na + ion within the structure, with f = 1 indicating a true random walk and f = 0 signifying truly non-random walk. Thus, f is different from HR, which describes the correlation between the motions of different Na + ions. Similarly, to our observations with HR, we find that correlation effects, as measured by f, are quite pronounced at all Na concentrations in NASICON (Figure 2d). For example, f decreases progressively from 0.2 to 0.01 as x varies from 1 to 2. Notably, f exhibits more significant non-monotonic behavior compared to HR in the range 0 < x < 3.
One of the main reasons behind the high degree of correlation (as indicated by HR and f) in NASICON is the high sensitivity of the E to the local Si/P environment (Figure 1e). Indeed, one may observe large variations in the local migration units, which leads to significant changes in E , at a given Na concentration. This is clearly shown in Figure S2, where the migration barrier varies between ~ 300 meV and ~ 800 meV at x ~ 2. Every Na + migration event leads to a dynamic change in local environment by creating new vacancies and eliminating existing ones, which can have sizeable effects on the migration barriers (and hence E ), eventually resulting in highly correlated motion.  (1) and Na (2), are shown as blue dashed lines. In panel b, the shaded regions mark areas where Na + -ion conductivity is highly (red) or slightly (blue) impacted by the charge carrier concentration.
It is worth analyzing how the distribution of Na-ions, across the Na(1) and Na(2) sites controls the Na diffusivity (and conductivity) of NASICON. Figure 3a shows the averaged Na occupancy of the Na(1) and Na(2) sites, as extracted from our kMC simulations. We compare these values with single-crystal X-ray diffraction experiments by Boilot et al. 30 at specific compositions, and our previous grandcanonical Monte Carlo simulations, 19 which do not account for the dynamics of the Na + ions.
In NASICONs, 18,48 Na-transport can only occur if multiple exchanges between Na(2)⟷ Na(1) sites occur. 18,30,31 This implies that both Na(1) and Na(2) must be occupied to some degree 18,30,31 for facile transport. The Na site occupancies derived from our kMC data which includes possible Na-migration events, at 443 K (170 °C, blue shapes in Figure 3a) show that both Na(1) and Na (2) are partially occupied at compositions 1 < x < 3. The kMC results are in qualitative agreement with X-ray measurements at 443 K of Ref. 30 (orange shapes in Figure 3a). In contrast, occupations generated by the grand-canonical Monte Carlo model (red filled shapes in Figure 3a) indicate the equilibrium Na occupancies, which are distinctly different from the kMC simulations (and experiments) that include the impact of Na movement.
Importantly, both our kMC and grand-canonical Monte Carlo simulations indicate that at Na1Zr2P3O12, Na + ions will be almost exclusively located at Na(1), even after including Na + dynamics at 443 K (Figure 3a). 32 In Na1Zr2P3O12, the Na(2) sites are empty because they are not thermodynamically stable at this composition. 32,41 Thus, the lack of Na atoms present in Na(2) sites that can act as diffusion carriers, contributes significantly to the low D (and σ) observed at low Na content (x ~ 0, Figure   2), in addition to the high E observed (i.e., at low values of Na@ Na(2) and Si@ Si/P in Figure 1e).
By contrast, at high Na contents (x ~ 3, Na4Zr2Si3O12), Na + diffusion is mostly limited by a lack of Na vacancies (either on Na(1) or Na(2) sites), which contributes to the moderately low computed D and σ. Note that the computed E barriers achieve their lowest values corresponding to high Na and Si content (Figure 1e), which is similar to the local migration units that will be observed at x ~ 3, and should, in principle, contribute to increase D and σ. However, the sharp decrease in the number of available vacancies prevails over the decrease in E , eventually lowering D and σ. While we only consider vacancy-based migration mechanisms in our work, a recent study has also explored interstitial-based migration mechanisms in various solid electrolytes. 49 Given the robust agreement between our predictions and measurements, we believe that an interstitial-based mechanism is not active within NASICON structures. followed by an eventual decline of both quantities as x → 3. In Figure 3a, at x = 3, both Na(1) and Na(2) sites are filled and no free charge carriers are available.
Experimentally and theoretically, [17][18][19]30 at x ~ 2.4 both Na(1) sites and Na (2) sites are expected to be occupied to some degree and the Na + has enough vacancies to migrate. Furthermore, we expect the overall migration barriers to be low at x ~ 2.4 compared to x < 2, as quantified by trends in E in Figure 1e. Therefore, Na3.4Zr2Si2.4P0.6O12 strikes the optimal combination of density of diffusion carriers (i.e., availability of Na ions and vacancies in both Na(1) and Na(2) sites), and low migration barriers to achieve the fastest Na-transport in the NASICON electrolyte.
We have synthesized different compositions of the NASICON, namely, Na2.5Zr2Si1.5P1.5O12, Na3Zr2Si2P1O12, and Na3.4Zr2Si2.4P0.6O12 and tested their Na-ion transport characteristics. Figure 4a shows the powder X-ray diffractions (PXRD) measured at room temperature, with the PXRD profiles matching the rhombohedral  The refined lattice parameters of Na2.5Zr2Si1.5P1.5O12, Na3Zr2Si2P1O12 and Na3.4Zr2Si2.4P0.6O12 are listed in Table 1 and are in agreement with existing reports. 17,50 The phase purity of these compounds enables us to characterize their variable temperature Na ionic conductivities, as reported in the Arrhenius plots of Figure 4b.
Impedance spectra at variable temperatures and their fitting details are also provided in Figure S9 and related text in the SI. The ionic conductivities of the synthesized Na1+xZr2SixP3-xO12 phases at variable temperature were fitted with appropriate equivalent circuits (see Methods section).  in the SI and includes the attempt frequency and the migration barrier. From our kMC data, we plot heatmaps indicating the spatial migration of Na-ions in Figure 5, along with the frequency associated with each Na migration event. Figure 5 shows the concentrations. Each panel in Figure 5 contains 1,024 migration units (Figure 1e), where each migration unit is represented by a distinct colored hexagonal prism. In  The direct visualization of the hopping frequency of Na + migration is insightful for unravelling the respective effects played by PO4 and SiO4 moieties on the overall Na + transport. In Figure 5, at low Na concentrations, e.g., x = 0.3, most migration units appear purple (low frequency) signifying low diffusivity/conductivity, whereas these blue regions turn progressively to yellow and eventually red for increasing values of x.
At x ~ 2, most Na-ions are actively engaged in ion transport as most migration units are red (Figure 5), highlighting facile ionic transport, consistent with the highest values of observed D and σ at x ~ 2.4 (Figure 2). As x approaches 3 (x ~ 2.97 shown in Figure 5), the color of the migration units eventually become purple, coinciding with the drop in D/σ at high Na-contents (Figure 2) due to the low availability of Navacancies. This informative analysis can be extended to other topical solid electrolytes with mixed polyanion sulfides, such as, GeS4, PS4, SiS4, etc. 8,9,[12][13][14] The underlying structural models to generate the data in Figure 5 are based on disordered NASICON structures, where at a given Na concertation the SiO4 and PO4 units are randomly distributed. This is also the case of the most researched sulfide electrolytes for Li-ion solid-state batteries with various degree of mixing of GeS4, SiS4, SiS4 and SnS4 moieties. 9,11-14 Thus, these NASICON structures truly mimic the experimental synthesis conditions (~ 1200 °C for solid-state synthesis 18,33 ) and/or heat treatments (e.g., low-temperature sol-gel synthesis 27,68 followed by sintering at ~1200 °C for densification 18,33,68 ), which quench the disorder of the phosphate and silicate units accessed at high temperatures. Previously we have shown that under equilibrium conditions NASICON should phase separate into P-rich and Si-rich domains, particularly across Na concentrations from x = 0 to 2. 19 It is also important to understand the impact of phase separation (or lack thereof) on Na-diffusivity and conductivity.
We have extended our kMC model to structures produced in regimes of complete phase separation in Figure S7 of the SI. Notably, the computed ionic conductivity and other related quantities predicted in regimes of phase separation look similar to those reported in Figure 2, with the exception of a dip in D and σ around x ~ 1.8 -2 compared to the scenario of fully disordered system (non-phase separated).
However, this dip in ion transport properties approximately accounts for one order of magnitude, well within the error of the available experiments and our calculations, thus making it challenging to detect the signature of phase separation via diffusivity/conductivity measurements or calculations. Our calculations suggest that Na diffusivity and conductivity cannot be used to identify the underlying phases within the NASICON electrolyte. We believe this conclusion may also be valid for other solid electrolytes and electrodes of interest, particularly those adopting the NASICON structure, which can also thermodynamically favor phase separation. 69,70 Notably, heatmaps produced in structures with phase separation (Figure S8) do clearly depict the presence of phase boundaries (differentiating Si-rich and P-rich domains), indicating that robust structural characterizations are still the best way of detecting phase separation.
Our computed barriers and kMC simulations suggest that higher Na + conductivity in NASICON can simply be achieved by increasing the SiO4 to PO4 ratio.
There are two reasons for this: (i) higher Si 4+ content decreases that of P 5+ cations in the system, thereby reducing the electrostatic repulsion between the Na + and P 5+

Sodium Migration Barriers from Density Functional Theory
Na migration barriers were calculated using the nudged elastic band (NEB) method 35 through density functional theory (DFT) simulations as implemented in the Vienna ab initio simulation package (VASP). 77,78 All the calculation parameters used in the DFT simulations can be found in Section S2 of the SI. The NEB barriers were modeled at 3 representative Na concentrations, at x = 0, 2, and 3 as in Na1+xZr2SixP3−xO12. This selection is motivated by our previous knowledge on the compositional phase diagram of Na1+xZr2SixP3−xO12, 19 where NASICON exhibits three distinct ground states at 0 K.
Special care is required for the NEB barrier calculations of NaZr2(PO4)3 (x = 0) which follows a pathway Na(1)→Na(2)→Na(1). In NaZr2(PO4)3 all the Na(1) sites are occupied and Na migration is only possible if vacancies in Na(1) are introduced. For this specific case the energy required for Na + migration should also include the formation energy of Na vacancy, which we computed using the method of Ref. 79 . We found that the Na vacancy formation energy in NaZr2(PO4)3 is ~ 474 meV, which we added in our model (see SI). Whereas the Na vacancy formation energy at x = 2 and 3 is negligible (~ 13 meV). Note that at x = 2, Na(1) and Na(2) are only partially occupied which guarantees facile Na transport.
To remove the directional dependence of migration barriers, we define the kinetically resolved activation (KRA) barrier, E of Equation 2, where the E barrier [Na(1) ⟷ Na (2)] is the NEB barrier, 39,40 and ∆Eend is the absolute difference between the computed energies of the initial and final end point structures.
The distribution of E at different configurations of Na and Si/P environments is shown in Figure S1, S2, S3, and Table S1 of the SI.

Local Cluster Expansion Hamiltonian
The calculated E are fitted to a local cluster expansion (LCE) Hamiltonian built around a migration unit (see Figure 1d) centered on the Na(1) site and using a cut-off radius of 5 Å. Details of the LCE model and the fitting strategy are given in the SI. The fitted Hamiltonian can be used efficiently to compute migration barriers inside the migration unit at any given Na/vacancy and Si/P content/configuration. The LCE Hamiltonian includes 1 point, 5 pair, and 1 triplet terms (see Table S3 and Table S4), respectively. The LCE Hamiltonian can reproduce the migration barriers with a RMS error of ~ ± 38 meV. The robustness of the LCE Hamiltonian was cross-validated using the leave-one out method. Notably, a variability of ~ ± 38 meV corresponds to a less than an order of magnitude in diffusivity (± 60 meV), 80 which is the typical uncertainty of experimental measurements and theoretical calculations for electrode materials. 40

Kinetic Monte Carlo Simulations
We implemented a rejection-free kinetic Monte Carlo (kMC) simulation scheme in an in-house-code, as described in Refs. 39,40 . The initial configurations for the kMC were generated using the canonical Monte Carlo simulation (based on our previous work where we constructed a global cluster expansion model on the NASICON) 19 in an 8×8×8 supercell at 0 ≤ x ≤ 3 containing 4,096 distinct Na sites and 3,072 Si/P sites.
These structures were generated at 1500 K, and close to the experimental synthesis temperature. 17 Variable temperature impedance measurements of Na1+xZr2SixP3−xO12 For the preparation of impedance measurements, the three samples were cold pressed into 6 mm pellets and sintered at ~ 1373 K for 12 hours. The sintered pellets were then gold-sputtered (~ 75 nm), and carbon paper (papyex) was added on each side of the pellet to ensure good contacts for the impedance measurements in the high temperature impedance cell BioLogic HTSH-100. The latter was inserted in high temperature furnace BioLogic HTF-1100 and connected to impedance analyzer The fitting of Nyquist plot was carried out using ZView software. For x = 1.5, the equivalent circuit used was based on semicircle phenomena as seen in other reports, 28,59,[84][85][86] with R1 representing the contact resistance, parallel of CPE1 and R2 representing the semi-circle, and initially a CPE2 for the diffusion tail. Since CPE2 has an n-value close to 0.5, a Warburg element is used to better fit the diffusion tail.

Conflicts of Interests
There are no conflicts to declare.

Supplementary information
The supplementary Information contains details of the density functional theory simulations, nudged elastic band results, kinetic Monte Carlo simulations, AC impedance measurements, equivalent circuits for the analysis of the AC-impedance data, and the Arrhenius plot of conductivities vs 1/temperatures of selected NASICON compositions.

Code and Data Availability
All the computational data associated with this study are freely available at the Git repository https://github.com/caneparesearch/NASICON_KMC_paper_data