A comprehensive exploration of Na⁺ ion transport in NaSICONs using molecular dynamics simulations

Abstract

Scandium-substituted Na_1+x+ySc_yZr_2−ySi_xP_3−xO₁₂ NaSICONs have emerged as promising electrolyte materials for all-solid-state sodium batteries. However, the comprehensive investigation of these multi-element structures is challenging due to their vast compositional space, leading to a limited number of compositions explored thus far. In this study, we address this issue by employing low-cost, yet high-precision force field molecular dynamics simulations based on density functional theory to investigate the Na⁺ mobility and resulting conductivity in Na_1+x+ySc_yZr_2−ySi_xP_3−xO₁₂ (0 ≤ x ≤ 3; 0 ≤ y ≤ 2). Our findings show that the incorporation of Sc³⁺- and Si⁴⁺-substituents enhances the conductivity, achieving values of 10⁻² S cm⁻¹ at room temperature for moderate to high substitution degrees. Moreover, our study demonstrates the efficacy of the applied methodology for large-scale screening, enabling the exploration of extensive configurational spaces of NaSICONs and other materials for potential use as solid-state electrolytes.

---

1. Introduction

Sodium-based energy storage systems have a longstanding history, but it is only in recent years that they stepped into the spotlight as concerns about limited resource availability of the market-dominating lithium-ion batteries prompted exploration of alternative solutions.^1–4 Sodium-ion batteries (NIBs) hold great promise as they provide a broad range of new possibilities, particularly for large-scale energy storage applications owing to the similar chemical properties of sodium and lithium, the abundance of sodium resources, and potential cost reduction.^1–15 Meanwhile, all-solid-state batteries (ASS) are being considered as the next generation of batteries due to their prolonged lifetimes, wide temperature ranges, and high energy and power densities.^16–26

The performance of ASS batteries depends crucially on the applied solid-state electrolyte (SSE), which must provide high conductivity at ambient temperatures. Among Na⁺-based SSEs, sodium superionic conductors (NaSICONs) have garnered significant attention, given their good thermal and (electro-)chemical stability, high ionic conductivity, and promising performance for ASS-NIBs.^{17–22,27–40} Several articles provide an overview of the developments of NaSICONs.^41–49 The parent NaSICON composition Na_1+xZr₂Si_xP_3−xO₁₂ (NZSP) exhibits the rhombohedral (R3̄c) structure above 150 °C, but undergoes a monoclinic (C2/c) phase transformation for 1.6 ≤ x ≤ 2.4 at room temperature.^50–58 As shown in Fig. 1a, the structure consists of a rigid three-dimensional framework composed of corner-sharing ZrO₆-octahedra and PO₄-tetrahedra, between which Na1 (6b) and Na2 (18e) sodium ion sites are located. In the non-substituted structure (x = 0), the energetically more favorable Na1 sites are fully occupied. As the P⁵⁺ sites are substituted by Si⁴⁺, the Na2 site occupation increases to satisfy the overall charge neutrality.^{50,51,54,55,57–67} The diffusion of the Na⁺ charge carriers is based on the correlated knock-off mechanism which involves the potential energy transfer from higher energy Na2 ions to lower energy Na1 ions.⁶⁸ This diffusion mechanism has already been proposed and studied in detail for various Na⁺ and Li⁺ conductors, such as NaSICONs,^68–70 garnets,⁶⁸ and Li₂CO₃.⁷¹ As illustrated in Fig. 1b, a Na2 ion jumps onto an adjacent occupied Na1 site and pushes the Na1 ion onto a neighboring unoccupied Na2 site. The bottlenecks of the migration are formed by opposing ZrO₆-octahedra and PO₄-tetrahedra adjacent to the pathway.^{50,51,55,60,61,63,65,72–79}

Fig. 1 (a) The primitive rhombohedral structure of NaSICONs with composition Na_1+xZr₂Si_xP_3−xO₁₂. (b) The correlated Na⁺ pushing-out conduction mechanism, indicated by arrows. One Na2 ion jumps onto an occupied Na1 site pushing the Na1 ion onto an adjacent unoccupied Na2 site. Key: Na1 (dark green), Na2 (light green), ZrO₆ (yellow), PO₄/SiO₄ (blue).

The three-dimensional Na⁺ mobility results in high ionic conductivities, e.g. up to 5.2 × 10⁻³ S cm⁻¹ for Na_3.4Zr₂Si_2.4P_0.6O₁₂ at room temperature.³³ However, NZSP is associated with a number of challenges. First, obtaining stoichiometric pure compositions is difficult due to the formation of secondary phases, such as ZrO₂. Second, the microstructure depends strongly on the sample preparation method. Both aspects have significant impact on the ionic conductivity.^21,45,80

Thus, numerous efforts in the past have been made to further improve NaSICONs in competition with liquid electrolytes. An outstanding feature of NaSICONs with the general formula NaM₂(AO₄)₃ is the ability of accommodating a wide range of M and A cations in the structural framework,⁴⁹ making multi-element substitution a feasible strategy for enhancing the performance of NaSICONs. For example, the substitution of Zr⁴⁺ by Sc³⁺ in NZSP leads to the suppression of secondary phases, optimization of the microstructure, and thus high Na⁺ conductivities.^81–94 Ma et al. reported a remarkable total conductivity of 4.0 × 10⁻³ S cm⁻¹ for Na_3.4Sc_0.4Zr_1.6Si₂P₁O₁₂ at room temperature.⁸¹ In addition, Sc³⁺-substituted structures exhibit a wide electrochemical stability window, suggesting both the application of high voltage cathode materials and the direct utilization of sodium metal anodes to boost the energy density.^35,81

However, the multi-element substitution of NaSICONs results in a vast compositional search space that requires fast, low-cost methods for investigating Na⁺ mobility and resulting conductivity. In this regard, computational techniques, such as molecular dynamics (MD) simulations, have proven as a powerful tool. MD provides comprehensive structural information and the dynamic evolution of conducting systems without relying on prior assumptions about the migration mechanism and pathways. Ionic jump events are directly observed from the modelled real-time dynamics and trajectories of the ions are used to identify the diffusion mechanism and quantify the diffusion properties. However, MD does require the accurate description of interatomic interactions, i.e. the potential energy surface (PES).^95,96

In ab initio MD (AIMD) simulations, the interactions between atoms are described by first-principles methods such as density functional theory (DFT), providing high reliability and accuracy as well as great transferability. However, due to its high computational expenses, AIMD is only feasible for small systems in terms of length scale (<1000 atoms) and timescale (∼100 ps) limiting the number of diffusion events.^95,96 As a result, diffusion properties estimated from AIMD simulations are often compromised by poor statistics and require high simulation temperatures to ensure a sufficient number of ion jumps.^95,97

In contrast, classical force field MD (FFMD) simulations employ interatomic potentials to model the interactions between the ions. The PES is described as a function of the atom positions rather than a variational ground state energy as in AIMD, which significantly reduces computational demands. Consequently, FFMD simulations are amenable for larger systems with up to millions of ions and longer simulation times (>1 ns), resulting in lower statistical variances in the determined diffusion properties. However, the accuracy of FFMD simulations depends on the availability and suitability of the FF parameter sets, which are usually restricted to a limited number of compounds and properties.^95,96

To tackle the trade-off between accuracy and efficiency in MD simulations, Kobayashi et al. introduced a high-throughput method for producing reliable FFs using a metaheuristic Cuckoo Search (CS) approach based on structure information from first-principles calculations. In this way, low-cost and accurate FFMD simulations are ensured.^98,99

It should be noted that data-driven machine learning potentials (MLP) has emerged as a powerful tool for materials exploration and development, providing an alternative to AIMD and classical FFMD simulations.^100–104 In particular, Neural Network Potentials (NNPs) have been demonstrated to reproduce AIMD results with high accuracy.^105,106 However, NNPs evaluate energy and forces in a purely empirical manner, which requires large amounts of training data. Additionally, their use in extrapolated spaces, beyond the training set, presents challenges. In contrast, the FF employed in this study, based on a physical model, is expected to maintain a higher degree of consistency even in extrapolated spaces. Moreover, NNPs face limitations in terms of computational time. For instance, in this study, the simulation of 1 ns (10⁶ steps) of a system with approximately 3500 atoms requires about 12 hours on 2 CPU cores. By comparison, performing the same calculation using M3GNet.py,^107,108 which is known for its relatively fast performance,^109,110 along with the Atomic Simulation Environment (ASE)¹¹¹ on 16 CPU cores, requires 458 hours. Thus, FFMD simulations are considered advantageous for studies that comprehensively evaluate a large number of compositions at multiple temperatures.

Here, DFT based FFMD simulations of Na_1+x+ySc_yZr_2−ySi_xP_3−xO₁₂ (0 ≤ x ≤ 3; 0 ≤ y ≤ 2), an example of a promising NaSICON compound, are performed to comprehensively investigate the Na⁺ transport behavior in multi-substituted NaSICONs. Initially, the FF parameters are optimized with the CS algorithm leveraging structural information acquired from short AIMD simulations of a training system. We show that the FF parameters obtained for the training system with a specific composition can be transferred to target systems with arbitrary compositions. The optimized potential parameters are then employed in the subsequent FFMD simulations across the entire compositional range of Na_1+x+ySc_yZr_2−ySi_xP_3−xO₁₂. Na⁺ diffusion coefficients are derived from the collective diffusion events described by their mean squared displacement. In addition, estimated radial distribution functions are analyzed to understand the influence of Sc³⁺- and Si⁴⁺-substitutions on the Na⁺ transport behavior.

On the one hand, this study aims to deepen the understanding of the influences of multi-element substitution on the transport of Na⁺ in NaSICONs, all in pursuit of achieving high ionic conductivity for ASS-NIBs applications. On the other hand, this study showcases the efficacy of the applied methodology for large-scale screening, enabling the exploration of extensive configurational spaces of NaSICONs and other materials promising as SSEs.

2. Methods

2.1 Force field model

FFMD simulations of Na⁺ diffusion in Na_1+x+ySc_yZr_2−ySi_xP_3−xO₁₂ require an FF model that describes the ion pair interactions and structural framework of NaSICONs. Thus, the FF model proposed by Kobayashi et al.⁹⁹ is applied, which consists of two-body potentials^112,113 and three-body angular potentials.¹¹⁴

The Morse potential is applied to describe the attractive interactions and short-range repulsions of cation–anion pairs,¹¹²

where D_0,ij is the equilibrium well depth defined by dissociation energy, α_ij describes the width of the potential associated to the bond softness, r_ij is the distance and r_0,ij the equilibrium distance between the i^th and j^th ion.

The screened Coulomb potential is added to account for the repulsive interactions among cations or anions. Since the Coulomb potential prevents cation or anion pairs with a small distance, their short-range repulsion can be neglected,¹¹²

where ε₀ is the permittivity of the vacuum, q_i1 and q_i2 are the formal charges of ions i₁ and i₂, and r_i1i2 is their ionic distance. The included complementary error function is a function of the ion distances over the screening length ρ_i1i2 = f(r_i1 + r_i2), which is determined by the sum of the effective radii r_i1 and r_i2 of the ions i₁ and i₂ and the screening factor f defining the radius of the Coulomb potential.

To avoid discontinuities in the potential energy and force, the following smoothing function with the cutoff radius r_c^2b is applied to both the Morse and Coulomb pair potential:

In addition, the Stillinger–Weber (SW) angular potential is applied to describe the framework of the non-migrating cations and oxygen anions,¹¹⁴

where r_ij and r_jk are the ionic distances and θ_ijk is the angle between the i^th, j^th, and k^th ion. λ_0,ijk is the equilibrium depth of the potential well and γ_0,ijk is the equilibrium angle between the i^th, j^th and k^th ion. The cutoff radius r_c^3b ensures that only nearest neighbor configurations are considered.

In previous studies, this FF model has been successfully applied to other solid state materials, such as LiZr₂(PO₄)₃,^99,115 La_(1−x)/3Li_xNbO₃,¹¹⁶ LiAl₅O₈,¹¹⁷ Li₃MCl₆,¹¹⁸ and La_1−xM_xF_3−x,⁹⁹ demonstrating rapid FF production with sufficient accuracy, as described below.

2.2 Force field parameterization

The FF parameters D_0,ij, α_ij and r_0,ij of the Morse potential, ρ_i1i2 of the Coulomb potential, and λ_ijk and γ_ijk of the SW potential must be determined for the Na_1+x+ySc_yZr_2−ySi_xP_3−xO₁₂ system. Thus, FF parametrization was conducted using a training system with 120 atoms according to the following steps:

(i) AIMD simulations were performed to obtain the equilibrium volume (V₀), the radial distribution function (RDF) and the bond angular distribution function (ADF) by averaging the ionic configurations of every 10^th step, omitting the first 800 steps to ensure that the system is in equilibrium.

The RDF g_ij(r) describes the particle density as a function of the distance r from the reference atom i to a surrounding atom j according to eqn (5),

with N_i and N_j being the numbers of the ions i and j, respectively.

The ADF h_ijk(θ) describes the particle density as a function of the bond angle θ_ijk between the reference atom i and two other atoms j and k within a given distance (2.8 Å) according to eqn (6),

with N_i, N_j and N_k being the numbers of the ion i, j and k, respectively.

Both RDF and ADF are determined by averaging the same element combinations, such as Sc–O, Zr–O, Si–O, P–O, and O–Sc–O, O–Zr–O, O–Si–O, O–P–O for Na_1+x+ySc_yZr_2−ySi_xP_3−xO₁₂. The RDF, ADF and V₀ are used as reference data in the subsequent FF fitting process.

(ii) The metaheuristic optimization of the FF parameters was performed by the CS algorithm^119–121 employing the AIMD reference data obtained in step (i).^98,99 The hyperparameters in the CS algorithm are the number of individuals (nests) in each generation (N_n) and the fraction of individuals replaced in the next generation by nests with new random solution (P_n).¹²¹

(iii) To evaluate the optimization process, FFMD simulations were performed using the obtained FF parameters. The results of the simulations were analyzed in the same way as for AIMD in step (i). The optimization process was checked for consistency using the loss function L assessing the differences between the RDF, ADF, and V₀ values obtained by FFMD and AIMD simulations (eqn (S1)–(S4) in ESI†).¹²²

Step (ii) and (iii) were repeated until the loss function value met the defined convergence criteria.^98,99

2.3 Force field molecular dynamics simulations and analysis

With the optimized FF parameters, FFMD simulations were performed for a total of 119 compositions of Na_1+x+yZr_2−ySc_ySi_xP_3−xO₁₂ (0 ≤ x ≤ 3; 0 ≤ y ≤ 2) in the temperature range of 298 ≤ T/K ≤ 798 in 100 K steps.

From the FFMD simulations, the mean squared displacement MSD of the migrating Na⁺ (eqn (7)) is obtained, which is used to determine the tracer diffusion coefficient D* (eqn (8)):

where r_i(t) and r_i(0) represent the position of ion i at time t and at t = 0, respectively, N is the number of mobile ions, d is the dimension of diffusion. D* describes the macroscopic ion transport and is related to the long-range diffusion coefficient D_σ by the Haven ratio H_R, which indicates the correlation of the trajectories of the migrating ions:^95,96

D* = H_RD_σ (9)

D_σ is described by the Nernst–Einstein relation:where u_i is the mobility and q_i is the charge of ion i, T is the temperature, and k_B is the Boltzmann constant. The mobility of the ion i characterizes its conductivity σ_i = u_iq_ic_i with c_i being the ion concentration.

Assuming H_R = 1, which is expected to have minimal impact on the results, the ionic conductivity is determined using the MSD, and the activation energy E_A is derived from the Arrhenius relationship:

Hereafter, D* will be referred to as D for simplicity.

2.4 Computational details

2.4.1 Ab initio molecular dynamics. NpT-AIMD simulations of the training system Na₃ScZrSiP₂O₁₂, containing 120 atoms, were performed for 7 ps with 1 fs time steps at 298 K to collect the data for the FF parameter fitting. The calculations were conducted based on Kohn–Sham DFT^123,124 using the Vienna Ab initio Simulation Package (VASP).^125–128 In the AIMD simulations the Parrinello–Rahman method with the Langevin thermostat was employed.^129,130 The generalized gradient approximation (GGA) functional parametrized by Perdew–Burke–Ernzerhof (PBEsol)¹³¹ was used. The projected augmented wave method (PAW)^132,133 was applied and plane waves¹²⁸ with an energy cut-off of 520 eV were employed. The electronic convergence criterion was set to 10⁻⁵ eV. To sample the Brillouin zone, a 2 × 2 × 2 gamma centered k-point mesh was used, following the Monkhorst and Pack scheme.¹³⁴ The 2p⁶3s¹ electrons of sodium, 2s²2p⁴ electrons of oxygen, 3s²3p⁶3d¹4s² electrons for scandium, 4s²4p⁶4d²5s² electrons of zirconium, 3s²3p² electrons of silicon and 3s²3p³ electrons of phosphorus were treated as valence electrons.

2.4.2 CS optimization. FF parameters of the Morse, Coulomb and SW potentials were determined using the Nagoya Atomistic Simulation Package (NAP).^98,99 The two- and three-body cutoff radii were set to r_c^2b = 7.0 Å and r_c^3b = 3.0 Å, respectively. The hyperparameters of the CS optimization were set to N_n = 36 and P_n = 0.25.

2.4.3 Force field molecular dynamics. All FFMD simulations were performed in cells with random Sc/Zr and Si/P arrangement using NAP.^98,99 Initial NpT-FFMD simulations of the training system Na₃ScZrSiP₂O₁₂, containing 120 atoms, were performed for 20 ps with 2 fs time steps at 298 K to validate the CS-derived FF parameters.

With the optimized FF parameters, FFMD simulations of supercells with composition Na_1+x+yZr_2−ySc_ySi_xP_3−xO₁₂ (0 ≤ x ≤ 3; 0 ≤ y ≤ 2) with up to 4013 atoms were performed in the temperature range of 298 ≤ T/K ≤ 798 in 100 K increments. First, NpT-FFMD simulations were conducted for 10 ps with 1 fs time steps to optimize the cell volume. Subsequently, NVT-FFMD simulations were performed for 1 ns with 1 fs time steps.

The Langevin thermostat, with a relaxation time of 50 fs, was used to control the temperature in both the NpT and NVT ensembles. For pressure regulation in the NpT ensemble, the Berendsen barostat with isotropic coupling and a relaxation of 50 fs was employed. To evaluate the influence of the thermostat on ionic diffusion, additional FFMD simulations were conducted under the NVE ensemble, omitting a heat bath. These simulations show that the diffusion coefficients increase by approximately a factor of 1.4 on average (Fig. S1†). Nevertheless, the overall trend of concentration-dependent diffusion behavior remains largely unchanged.

3. Results

3.1 Validation of optimized force field parameters

The optimization of the FF parameters of the six-component NaSICON structure Na_1+x+ySc_yZr_2−ySi_xP_3−xO₁₂ follows a multiple-step process detailed in Section 2 of the ESI.† The FF parameters were optimized using the training system Na₃ScZrSiP₂O₁₂. For further validation, the obtained FF parameters are applied in FFMD simulations across the entire compositional range of Na_1+xZr₂Si_xP_3−xO₁₂ (0 ≤ x ≤ 3) at 598 K and the resulting conductivities are compared with literature values. These compositions and temperature are chosen due to the extensive data available in the literature, allowing a comprehensive comparison. Fig. 2 shows the comparison of the resulting conductivities with literature data from experiments and KMC simulations.^{51,52,55,60,66,135} Although, the data in the literature exhibit significant variation that exceed the uncertainty of the simulation, the conductivities obtained from FFMD simulations in this work and the literature data show similar dependence on the substitution degree x. This indicates that the FF parameters not only adequately describe the NaSICON structure but also capture the interactions between the migrating Na⁺ and both the structural framework and the neighboring Na⁺, which are crucial for the Na⁺ transport in NaSICONs, as highlighted in our previous studies.^135,136 It is shown that the parameters originally obtained for the training system with the specific composition Na₃ScZrSiP₂O₁₂ can be applied to the entire compositional range of Na_1+x+ySc_yZr_2−ySi_xP_3−xO₁₂.

Fig. 2 Comparison of Na⁺ conductivity of Na_1+xZr₂Si_xP_3−xO₁₂ (0 ≤ x ≤ 3) at 598 K obtained by FFMD simulations using FF parameters derived from the CS optimization in the Na₃ScZrSiP₂O₁₂ training system with literature data obtained from experimental studies^51,52,55,66 and KMC simulations^60,135 at 573 K.

3.2 Microscopic Na⁺ diffusion

Since the optimized FF parameter adequately describe the NaSICON system, NVT-MD simulations are conducted using this FF parameter set (Table S3†) to investigate the Na⁺ diffusion behavior in Na_1+x+yZr_2−ySc_ySi_xP_3−xO₁₂ (0 ≤ x ≤ 3; 0 ≤ y ≤ 2) over the temperatures range of 298 ≤ T/K ≤ 798. Information about the microscopic Na⁺ diffusion behavior can be directly obtained from the MD simulations.

Fig. 3 depicts the resulting density plot of the accumulated trajectories of the Na⁺ in Na₃Sc_1.5Zr_0.5Si_0.3P_2.7O₁₂ over the simulation time, revealing the well-connected three-dimensional Na⁺ diffusion pathways between the Na1 and Na2 sites that enable the long-range Na⁺ transport. This conduction pathway is in excellent agreement with the findings of Deng et al.,¹³⁷ who used a combination of FFMD simulations with parameters of the interatomic potentials developed by Pedone et al.,¹³⁸ bond valence energy landscape method, maximum-entropy method/Rietveld analysis method.

Fig. 3 The density plot of the accumulated trajectories of Na⁺ ions in Na₃Sc_1.5Zr_0.5Si_0.3P_2.7O₁₂ over the MD simulation time. The close-up view shows the well-connected diffusion pathways between the Na1 (light green) and the Na2 (dark green) sites.

Fig. 4a depicts the resulting MSDs of the Na⁺ in Na₃Sc_1.5Zr_0.5Si_0.3P_2.7O₁₂ at all simulation temperatures. In Fig. 4b, the MSDs of all other ions in the same compound are shown for 698 K. As anticipated, solely the Na⁺ diffuse in the cell, whereas O²⁻, Sc³⁺, Zr⁴⁺, Si⁴⁺, P⁵⁺ vibrate around their initial equilibrium sites.

Fig. 4 MSD of (a) Na⁺ in the temperature range of 298 ≤ T/K ≤ 798 and (b) O²⁻, Sc³⁺, Zr⁴⁺, Si⁴⁺, P⁵⁺ at 698 K in Na₃Sc_1.5Zr_0.5Si_0.5P_2.5O₁₂.

3.3 Exploration of the configurational space

The microscopic Na⁺ diffusion behavior is accurately captured using the optimized FF parameter set. Consequently, the diffusion coefficient, conductivity, and activation energy are predicted across the entire configurational space of Na_1+x+ySc_yZr_2−ySi_xP_3−xO₁₂ (0 ≤ x ≤ 3; 0 ≤ y ≤ 2). As most current studies focus on only a few compositions, this comprehensive prediction offers a broader picture of Na⁺ diffusion in the Sc³⁺- and Si⁴⁺-substituted NaSICON system.

The Na⁺ tracer diffusion coefficients D are calculated using the MSDs according to eqn (8). Following the approach proposed by Usler et al.,¹³⁹ the error u(D) in D is estimated, indicating that u(D) is two orders of magnitude smaller than D for temperatures above 498 K. Fig. S9† illustrates the Arrhenius plots for various compositions with different substitution ratios and Na⁺ concentrations, showing the anticipated linear relationship between D and temperatures above 498 K. Below this temperature, however, deviations from the linearity are observed, probably due to insufficient statistics of Na⁺ hopping events. The Na⁺ conductivities σ and activation energies E_A are derived from eqn (11). For the determination of E_A, only diffusion coefficients obtained at temperatures between 498 ≤ T/K ≤ 798 are considered to ensure a statistically significant number of Na⁺ jumps. The calculated values of D and σ at 598 K, E_A and the cell volumes V of Na_1+x+ySc_yZr_2−ySi_xP_3−xO₁₂ (0 ≤ x ≤ 3; 0 ≤ y ≤ 2) are depicted in Fig. 5. Note that the Na⁺ diffusion coefficients and conductivities are shown for 598 K rather than room temperature because the conductivity values below 498 K fall within the extrapolated range. Additionally, 598 K is particularly relevant for practical applications, as it corresponds to the operating temperature of NAS batteries.¹⁴⁰

Fig. 5 (a) Na⁺ tracer diffusion coefficient D at 598 K, (b) Na⁺ conductivity σ obtained at 598 K, (c) activation energy E_A obtained from D at temperatures between 498 ≤ T/K ≤ 798, and (d) cell volume V as a function of the Sc³⁺-content y and Si⁴⁺-content x in Na_1+x+yZr_2−ySc_ySi_xP_3−xO₁₂.

It is noted that in this study, the Sc/Zr and Si/P arrangements were randomly configured, and thus the impact of such cation arrangements on the diffusion coefficient of Na⁺ are investigated. For ten randomly selected compositions, four random initial arrangements were created and the relative deviation of the diffusion coefficient among these was evaluated. The distribution of the obtained diffusion coefficients is shown in Fig. S10.† The variation in diffusion coefficients due to differences in arrangement are very small, with the relative standard error being a maximum of 5.6% and an average of 2.6% across the 10 compositions. Therefore, it is concluded that the supercell size and simulation time (1 ns) applied in this study are sufficiently large.

In our previous studies on Na_1+xM₂Si_xP_3−xO₁₂, we have revealed the dependence of the Na⁺ mobility in NaSICONs on both the local ionic environment of the Na⁺ lattice sites and the Na⁺ migration energy.^135,136 The introduction of Si⁴⁺-substituents on P⁵⁺ sites corresponds to a charged defect of −1, leading to favorable Na⁺–Si⁴⁺ pair interactions which influence the Na⁺ mobility. Three compositional ranges in Na_1+xM₂Si_xP_3−xO₁₂ can be distinguished:

(i) At small substitution concentrations (x < 0.5), there are only a few Na⁺ charge carriers, that are trapped in low energy states near the introduced Si⁴⁺ due to the attractive Na⁺–Si⁴⁺ interactions. The depth and number of traps increases with the number of Si⁴⁺, resulting in low conductivities and an increase of the activation energies.

(ii) At intermediate to high substitution concentrations (0.5 ≤ x ≤ 2.0), Na⁺ percolation pathways form between favorable Si⁴⁺-rich positions. The Na⁺ mobility is thus determined by the migration energy. As a result, the conductivity increases, and the activation energy decreases with substitution content.

(iii) At high substitution concentrations (x > 2.0), the insufficient number of vacant Na⁺ lattice sites lead to a decrease in ionic conductivity while the activation energy is the lowest due to the low migration energy. However, strong unfavorable Na⁺–Na⁺ repulsion may lead to an increase in activation energy.

Accordingly, the behavior of D (Fig. 5a), σ (Fig. 5b), and E_A (Fig. 5c) can be described as a function of the substitution degree in Na_1+x+ySc_yZr_2−ySi_xP_3−xO₁₂. In this system, the substitution of Sc³⁺ for Zr⁴⁺ and Si⁴⁺ for P⁵⁺ creates charged defects of −1. At minor substituents concentrations, both D and σ remain low, while E_A increases, although less pronounced than previously observed in Na_1+xM₂Si_xP_3−xO₁₂.¹³⁵ This behavior can be attributed to the limited number of Na⁺ charge carriers being trapped in low energy states near the introduced Sc³⁺ and Si⁴⁺ ions. With increasing substitution level, percolation pathways form between Sc³⁺- and Si⁴⁺-rich environments, leading to an increase in both D and σ and a decrease in E_A. At high substitution concentrations, D and σ decrease due to the lack of available vacant Na⁺ sites. But the activation energy is lowest, particularly in the single-substituted structures, as previously observed in Na_1+xZr₂Si_xP_3−xO₁₂. Though, at high Sc³⁺- and intermediate Si⁴⁺-content, the activation energy is elevated. This can be attributed to strong Na⁺–Na⁺ repulsive interactions that prevent the occupation of Na⁺ sites, as found in Na_1+xSn₂Si_xP_3−xO₁₂.¹³⁵ In summary, the simulations reveal that the region of high diffusivity and conductivity follows a diagonal pattern, with increasing values observed at high Sc content (1.3 ≤ y ≤ 1.9) and medium to high Si content (1.3 ≤ x ≤ 2.7).

The comparison of σ (Fig. 5b) and E_A (Fig. 5c) with cell volume V (Fig. 5d) shows that there is no clear relation, suggesting that the enlargement of the bottleneck area with increased cell volume has a less pronounced impact on conductivity and activation energy than previously proposed in the literature.^{49,50,55,65,74–76,78,79}

Earlier studies have already experimentally investigated some compositions of Na_1+x+ySc_yZr_2−ySi_xP_3−xO₁₂. For example, Lunghammer et al. employed time-domain nuclear magnetic resonance (NMR) and broadband conductivity spectroscopy to investigate the Na⁺ motion in the Si⁴⁺-rich composition Na_3.4Sc_0.4Zr_1.6Si₂PO₁₂ on both short- and long-range scales revealing activation energies of 0.13 eV to 0.15 eV and 0.29 eV to 0.31 eV, respectively.⁸² Rapid Na⁺ exchange leads to ionic conductivities of 2 × 10⁻³ S cm⁻¹ at 25 °C.⁸² For the same composition, Ma et al. reported bulk conductivities of 6.2 × 10⁻³ S cm⁻¹ at 298 K.⁸¹ For the Sc³⁺-rich composition Na_3.4Sc₂Si_0.4P_2.6O₁₂, Kaus et al. observed values of E_A = 0.2 eV and σ = 4 × 10⁻³ S cm⁻¹ at 350 K using NMR relaxometry.⁸³ In a separate study, Guin et al. conducted impedance measurements on the same composition, reporting an activation energy of 0.14 eV, which corresponds to the bulk transport.⁸⁴ In comparison, MD simulations of this study yield values of E_A = 0.11 eV and σ = 2.6 × 10⁻² S cm⁻¹ at 298 K for the Si⁴⁺-rich composition Na_3.4Sc_0.4Zr_1.6Si₂PO₁₂, and E_A = 0.08 eV and σ = 1 × 10⁻¹ S cm⁻¹ at 350 K for the Sc³⁺-rich composition Na_3.4Sc₂Si_0.4P_2.6O₁₂. The values obtained for E_A in this study are consistent with the values for short-range bulk Na⁺ transport, capturing the elementary diffusion processes involving local, correlated motions influenced by Coulomb interactions and structural disorder.⁸² However, the obtained values for σ are one to two orders of magnitude higher than those reported in previously cited studies.^81–84 Nevertheless, our findings are in good agreement with the conductivities determined by Subramanian et al. at 573 K, using impedance measurements for various compositions, such as σ = 1.35 × 10⁻¹ S cm⁻¹ versus σ = 1.07 × 10⁻¹ S cm⁻¹ for Na₃Sc_0.5Zr_1.5Si_1.5P_1.5O₁₂.⁹⁰ Moreover, similar general behavior of conductivity and activation energy as a function of the substitution concentration have been reported in the literature for Na_3+ySc_yZr_2−ySi₂PO₁₂ (0 ≤ y ≤ 0.6) and of Na_3+xSc₂Si_xP_3−xO₁₂ (0.05 ≤ x ≤ 0.8).^81,84,94

The effect of Sc³⁺- and Si⁴⁺-substitution on the Na⁺ diffusion behavior in NaSICONs is further clarified by estimating RDFs of cations adjacent to the Na⁺ sites. Fig. 6a illustrates der RDFs of Si⁴⁺–Na⁺ and P⁵⁺-Na⁺ in Na_1.5Zr₂Si_0.5P_2.5O₁₂. At smaller distances, the peak of the RDF of Si⁴⁺–Na⁺ pairs is larger compared to that of P⁵⁺–Na⁺ pairs, which indicates that more Na⁺ are adjacent to Si⁴⁺ than to P⁵⁺, suggesting the Na⁺ trapping by the Si⁴⁺-substituents.

Fig. 6 Comparison of (a) RDFs of Si⁴⁺–Na⁺ (solid line) and P⁵⁺–Na⁺ (dashed line) in Na_1.5Zr₂Si_0.5P_2.5O₁₂ and (b) RDFs of Sc³⁺–Na⁺ (solid line) and Zr⁴⁺–Na⁺ (dash line) in Na_1.3Sc_0.3Zr_1.7P₃O₁₂.

Similarly, Fig. 6b shows the comparison of the RDFs of Sc³⁺–Na⁺ and Zr⁴⁺–Na⁺ pairs in Na_1.3Sc_0.3Zr_1.7P₃O₁₂. The RDF of Sc³⁺–Na⁺ pairs exhibit larger peak than that of Zr⁴⁺–Na⁺ pairs at small distances, indicating that the Na⁺ are trapped by Sc³⁺ substituents.

4. Conclusion

The Na⁺ transport properties in NaSICONs were exhaustively investigated in Na_1+x+ySc_yZr_2−ySi_xP_3−xO₁₂ (0 ≤ x ≤ 3; 0 ≤ y ≤ 2) using FFMD simulations. For this purpose, FF parameters were determined based on the structural information of the training system Na₃ScZrSiP₂O₁₂ obtained from short AIMD simulations. The obtained FF parameters can be applied to target systems with arbitrary compositions, which they describe with performance comparable to DFT calculations. This underscores the effectiveness of the FF model in characterizing the structure and interactions in NaSICONs and in describing the Na⁺ transport behavior without necessitating explicit consideration of the ion migration and transition state structure. In this way, the exhaustive exploration of the Na⁺ transport in NaSICONs across a broad spectrum of compositions is feasible using fast, low-cost, and yet accurate FFMD simulations.

NVT-MD simulations revealed the three-dimensional diffusion pathways of Na⁺ ions between the Na1 and Na2 sites. The resulting Na⁺ diffusion coefficients D, ionic conductivities σ, and activation energies E_A show similar dependence on the substitution concentration as previously observed for the single-substituted structures Na_1+xM₂Si_xP_3−xO₁₂. The Na⁺ transport behavior is primarily governed by the effects of Na⁺ trapping and percolation caused by the introduced substituents. At minor substitution concentrations, Na⁺ are trapped in low energy states near the substituents, resulting in low values of D and σ and high values of E_A. At intermediate substitution levels, Na⁺ percolation pathways are formed, leading to an increase in D and σ and decrease in E_A. However, at very high substitution levels, D and σ decrease due to the lack of available Na⁺ sites. Overall, the simulations indicate that the highest diffusivity and conductivity occur along a diagonal compositional trend, corresponding to elevated Sc contents (1.3 ≤ y ≤ 1.9) and moderate to high Si contents (1.3 ≤ x ≤ 2.7). The activation energies are found to be in good agreement with experimental values describing short-range Na⁺ ion motions.

In addition, the Na⁺ trapping effects of the introduced Sc³⁺- and Si⁴⁺-substitution was verified by estimating RDFs, providing further insights into the Na⁺ migration behavior.

This study demonstrates the efficacy of the applied methodology in exploring large configurational spaces of NaSICONs and predicting promising materials for use in ASS-NIBs. The FFMD simulations, applying FF parameters optimized for the particular system under investigation, offer a computationally efficient way to identify substitution levels that lead to the highest conductivities. This allows further investigations to focus on this optimal range instead of exploring large compositional spaces. The key advantage of the applied method is its simplicity and versatility, making it applicable to a wide range of compounds and structures. In addition, the crucial impact of substitution in NaSICONs is further elucidated, indicating the high potential of these materials, especially Na_1+x+ySc_yZr_2−ySi_xP_3−xO₁₂.

Data availability

The data supporting this article have been included as part of the ESI.† The raw data is available from the authors upon request.

Author contributions

M. N., J. S. and S. N. conceived the study. J. S. and S. T. performed the simulations. All the authors contributed to the discussion and the paper writing.

Conflicts of interest

The authors declare no competing financial or non-financial interests.

Acknowledgements

This research was partially funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) – project number 452855747. J. S. is International Research Fellow of Japan Society for the Promotion of Science.

References

1. C. Delmas, Sodium and sodium-ion batteries: 50 years of research, Adv. Energy Mater., 2018, 8, 1703137.
2. N. Yabuuchi, K. Kubota, M. Dahbi and S. Komaba, Research development on sodium-ion batteries, Chem. Rev., 2014, 114, 11636–11682.
3. D. Larcher and J.-M. Tarascon, Towards greener and more sustainable batteries for electrical energy storage, Nat. Chem., 2015, 7, 19–29.
4. D. Kundu, E. Talaie, V. Duffort and L. F. Nazar, The emerging chemistry of sodium ion batteries for electrochemical energy storage, Angew. Chem., Int. Ed., 2015, 54, 3431–3448.
5. B. L. Ellis and L. F. Nazar, Sodium and sodium-ion energy storage batteries, Curr. Opin. Solid State Mater. Sci., 2012, 16, 168–177.
6. V. Palomares, P. Serras, I. Villaluenga, K. B. Hueso, J. Carretero-González and T. Rojo, Na-ion batteries, recent advances and present challenges to become low cost energy storage systems, Energy Environ. Sci., 2012, 5, 5884–5901.
7. M. D. Slater, D. Kim, E. Lee and C. S. Johnson, Sodium-ion batteries, Adv. Funct. Mater., 2013, 23, 947–958.
8. H. Pan, Y.-S. Hu and L. Chen, Room-temperature stationary sodium-ion batteries for large-scale electric energy storage, Energy Environ. Sci., 2013, 6, 2338–2360.
9. S. Y. Hong, Y. Kim, Y. Park, A. Choi, N.-S. Choi and K. T. Lee, Charge carriers in rechargeable batteries: Na ions vs. Li ions, Energy Environ. Sci., 2013, 6, 2067–2081.
10. J.-Y. Hwang, S.-T. Myung and Y.-K. Sun, Sodium-ion batteries: present and future, Chem. Soc. Rev., 2017, 46, 3529–3614.
11. A. Eftekhari and D.-W. Kim, Sodium-ion batteries: new opportunities beyond energy storage by lithium, J. Power Sources, 2018, 395, 336–348.
12. C. Vaalma, D. Buchholz, M. Weil and S. Passerini, A cost and resource analysis of sodium-ion batteries, Nat. Rev. Mater., 2018, 3, 1–11.
13. L. Chen, M. Fiore, J. E. Wang, R. Ruffo, D.-K. Kim and G. Longoni, Readiness level of sodium-ion battery technology: a materials review, Adv. Sustainable Syst., 2018, 2, 1700153.
14. T. Liu, Y. Zhang, Z. Jiang, X. Zeng, J. Ji, Z. Li, X. Gao, M. Sun, Z. Lin and M. Ling, Exploring competitive features of stationary sodium ion batteries for electrochemical energy storage, Energy Environ. Sci., 2019, 12, 1512–1533.
15. Y. E. Durmus, H. Zhang, F. Baakes, G. Desmaizieres, H. Hayun, L. Yang, M. Kolek, V. Küpers, J. Janek and D. Mandler, Side by side battery technologies with lithium-ion based batteries, Adv. Energy Mater., 2020, 10, 2000089.
16. C. Jiang, H. Li and C. Wang, Recent progress in solid-state electrolytes for alkali-ion batteries, Sci. Bull., 2017, 62, 1473–1490.
17. W. Hou, X. Guo, X. Shen, K. Amine, H. Yu and J. Lu, Solid electrolytes and interfaces in all-solid-state sodium batteries: Progress and perspective, Nano Energy, 2018, 52, 279–291.
18. C. Zhao, L. Liu, X. Qi, Y. Lu, F. Wu, J. Zhao, Y. Yu, Y.-S. Hu and L. Chen, Solid-state sodium batteries, Adv. Energy Mater., 2018, 8, 1703012.
19. C. Zhou, S. Bag and V. Thangadurai, Engineering materials for progressive all-solid-state Na batteries, ACS Energy Lett., 2018, 3, 2181–2198.
20. L. Fan, S. Wei, S. Li, Q. Li and Y. Lu, Recent progress of the solid-state electrolytes for high-energy metal-based batteries, Adv. Energy Mater., 2018, 8, 1702657.
21. Q. Ma and F. Tietz, Solid-state electrolyte materials for sodium batteries: towards practical applications, ChemElectroChem, 2020, 7, 2693–2713.
22. V. Thangadurai and B. Chen, Solid Li-and Na-Ion Electrolytes for Next Generation Rechargeable Batteries, Chem. Mater., 2022, 34, 6637–6658.
23. H. Che, S. Chen, Y. Xie, H. Wang, K. Amine, X.-Z. Liao and Z.-F. Ma, Electrolyte design strategies and research progress for room-temperature sodium-ion batteries, Energy Environ. Sci., 2017, 10, 1075–1101.
24. T. Famprikis, P. Canepa, J. A. Dawson, M. S. Islam and C. Masquelier, Fundamentals of inorganic solid-state electrolytes for batteries, Nat. Mater., 2019, 18, 1278–1291.
25. D. H. S. Tan, A. Banerjee, Z. Chen and Y. S. Meng, From nanoscale interface characterization to sustainable energy storage using all-solid-state batteries, Nat. Nanotechnol., 2020, 15, 170–180.
26. Y. Wang, S. Song, C. Xu, N. Hu, J. Molenda and L. Lu, Development of solid-state electrolytes for sodium-ion battery–a short review, NMS, 2019, 1, 91–100.
27. Y. Noguchi, E. Kobayashi, L. S. Plashnitsa, S. Okada and J. Yamaki, Fabrication and performances of all solid-state symmetric sodium battery based on NASICON-related compounds, Electrochim. Acta, 2013, 101, 59–65.
28. F. Lalère, J.-B. Leriche, M. Courty, S. Boulineau, V. Viallet, C. Masquelier and V. Seznec, An all-solid state NASICON sodium battery operating at 200 °C, J. Power Sources, 2014, 247, 975–980.
29. S. Song, H. M. Duong, A. M. Korsunsky, N. Hu and L. Lu, A Na⁺ superionic conductor for room-temperature sodium batteries, Sci. Rep., 2016, 6, 1–10.
30. W. Zhou, Y. Li, S. Xin and J. B. Goodenough, Rechargeable sodium all-solid-state battery, ACS Cent. Sci., 2017, 3, 52–57.
31. Z. Zhang, Q. Zhang, J. Shi, Y. S. Chu, X. Yu, K. Xu, M. Ge, H. Yan, W. Li and L. Gu, A self-forming composite electrolyte for solid-state sodium battery with ultralong cycle life, Adv. Energy Mater., 2017, 7, 1601196.
32. T. Lan, C.-L. Tsai, F. Tietz, X.-K. Wei, M. Heggen, R. E. Dunin-Borkowski, R. Wang, Y. Xiao, Q. Ma and O. Guillon, Room-temperature all-solid-state sodium batteries with robust ceramic interface between rigid electrolyte and electrode materials, Nano Energy, 2019, 65, 104040.
33. Q. Ma, C.-L. Tsai, X.-K. Wei, M. Heggen, F. Tietz and J. T. S. Irvine, Room temperature demonstration of a sodium superionic conductor with grain conductivity in excess of 0.01 S cm⁻¹ and its primary applications in symmetric battery cells, J. Mater. Chem., 2019, 7, 7766–7776.
34. P. Kehne, C. Guhl, Q. Ma, F. Tietz, L. Alff, R. Hausbrand and P. Komissinskiy, Electrochemical performance of all-solid-state sodium-ion model cells with crystalline Na_xCoO₂ thin-film cathodes, J. Electrochem. Soc., 2019, 166, A5328.
35. P. Kehne, C. Guhl, Q. Ma, F. Tietz, L. Alff, R. Hausbrand and P. Komissinskiy, Sc-substituted Nasicon solid electrolyte for an all-solid-state Na_xCoO₂/Nasicon/Na sodium model battery with stable electrochemical performance, J. Power Sources, 2019, 409, 86–93.
36. Y. Lu, J. A. Alonso, Q. Yi, L. Lu, Z. L. Wang and C. Sun, A high-performance monolithic solid-state sodium battery with Ca²⁺ doped Na₃Zr₂Si₂PO₁₂ electrolyte, Adv. Energy Mater., 2019, 9, 1901205.
37. E. Matios, H. Wang, C. Wang, X. Hu, X. Lu, J. Luo and W. Li, Graphene regulated ceramic electrolyte for solid-state sodium metal battery with superior electrochemical stability, ACS Appl. Mater. Interfaces, 2019, 11, 5064–5072.
38. Z. Zhang, S. Wenzel, Y. Zhu, J. Sann, L. Shen, J. Yang, X. Yao, Y.-S. Hu, C. Wolverton and H. Li, Na₃Zr₂Si₂PO₁₂: a stable Na⁺-ion solid electrolyte for solid-state batteries, ACS Appl. Energy Mater., 2020, 3, 7427–7437.
39. J. Yang, G. Liu, M. Avdeev, H. Wan, F. Han, L. Shen, Z. Zou, S. Shi, Y.-S. Hu and C. Wang, Ultrastable all-solid-state sodium rechargeable batteries, ACS Energy Lett., 2020, 5, 2835–2841.
40. X. Wang, W. Mei, J. Chen, D. Wang and Z. Mao, Rare Earth Oxide-Assisted Sintered NASICON Electrolyte Composed of a Phosphate Grain Boundary Phase with Low Electronic Conductivity, ACS Appl. Energy Mater., 2021, 5, 777–783.
41. J. Alamo and R. Roy, Crystal chemistry of the NaZr₂(PO₄)₃, NZP or CTP, structure family, J. Mater. Sci., 1986, 21, 444–450.
42. N. a. Anantharamulu, K. Koteswara Rao, G. Rambabu, B. Vijaya Kumar, V. Radha and M. Vithal, A wide-ranging review on Nasicon type materials, J. Mater. Sci., 2011, 46, 2821–2837.
43. R. Rajagopalan, Z. Zhang, Y. Tang, C. Jia, X. Ji and H. Wang, Understanding crystal structures, ion diffusion mechanisms and sodium storage behaviors of NASICON materials, Energy Storage Mater., 2021, 34, 171–193.
44. Z. Yang, B. Tang, Z. Xie and Z. Zhou, NASICON-type Na₃Zr₂Si₂PO₁₂ solid-state electrolytes for sodium batteries, ChemElectroChem, 2021, 8, 1035–1047.
45. Y. B. Rao, K. K. Bharathi and L. N. Patro, Review on the synthesis and doping strategies in enhancing the Na ion conductivity of Na₃Zr₂Si₂PO₁₂ (NASICON) based solid electrolytes, Solid State Ionics, 2021, 366, 115671.
46. C. Li, R. Li, K. Liu, R. Si, Z. Zhang and Y.-S. Hu, NaSICON: a promising solid electrolyte for solid-state sodium batteries, Interdiscip. Mater., 2022, 1, 396–416.
47. L. Zhang, Y. Liu, Y. You, A. Vinu and L. Mai, NASICONs-type solid-state electrolytes: the history, physicochemical properties, and challenges, Interdiscip. Mater., 2023, 2, 91–110.
48. K. Singh, A. Chakraborty, R. Thirupathi and S. Omar, Recent advances in NASICON-type oxide electrolytes for solid-state sodium-ion rechargeable batteries, Ionics, 2022, 28, 5289–5319.
49. M. Guin and F. Tietz, Survey of the transport properties of sodium superionic conductor materials for use in sodium batteries, J. Power Sources, 2015, 273, 1056–1064.
50. H.-P. Hong, Crystal structures and crystal chemistry in the system Na_1+xZr₂Si_xP_3−xO₁₂, Mater. Res. Bull., 1976, 11, 173–182.
51. J. B. Goodenough, H.-P. Hong and J. A. Kafalas, Fast Na⁺-ion transport in skeleton structures, Mater. Res. Bull., 1976, 11, 203–220.
52. J. P. Boilot, J. P. Salanie, G. Desplanches and D. Le Potier, Phase transformation in Na_1+xSi_xZr₂P_3−xO₁₂ compounds, Mater. Res. Bull., 1979, 14, 1469–1477.
53. U. von Alpen, M. F. Bell and W. Wichelhaus, Phase transition in nasicon (Na₃Zr₂Si₂PO₁₂), Mater. Res. Bull., 1979, 14, 1317–1322.
54. J.-J. Didisheim, E. Prince and B. J. Wuensch, Neutron Rietveld analysis of structural changes in NASICON solid solutions Na_1+xZr₂Si_xP_3−xO₁₂ at elevated temperatures: x = 1.6 and 2.0 at 320 °C, Solid State Ionics, 1986, 18, 944–958.
55. W. H. Baur, J. R. Dygas, D. H. Whitmore and J. Faber, Neutron powder diffraction study and ionic conductivity of Na₂Zr₂SiP₂O₁₂ and Na₃Zr₂Si₂PO₁₂, Solid State Ionics, 1986, 18, 935–943.
56. T. Feist, P. K. Davies and E. Vogel, The energetics of phase transitions in the system Na_1+xSi_xZr₂P_3−xO₁₂ 1.9 ≤ x ≤ 2.5, Thermochim. Acta, 1986, 106, 57–61.
57. J. P. Boilot, G. Collin and P. Colomban, Crystal structure of the true NASICON: Na₃Zr₂Si₂PO₁₂, Mater. Res. Bull., 1987, 22, 669–676.
58. J.-P. Boilot, G. Collin and P. Colomban, Relation structure-fast ion conduction in the NASICON solid solution, J. Solid State Chem., 1988, 73, 160–171.
59. L.-O. Hagman, P. Kierkegaard, P. Karvonen, A. I. Virtanen and J. Paasivirta, The crystal structure of NaMe₂^IV(PO₄)₃; Me^IV = Ge, Ti, Zr, Acta Chem. Scand., 1968, 22, 1822–1832.
60. Z. Deng, T. P. Mishra, E. Mahayoni, Q. Ma, A. J. K. Tieu, O. Guillon, J.-N. Chotard, V. Seznec, A. K. Cheetham, C. Masquelier, G. S. Gautam and P. Canepa, Fundamental investigations on the sodium-ion transport properties of mixed polyanion solid-state battery electrolytes, Nat. Commun., 2022, 13, 4470.
61. P. P. Kumar and S. Yashonath, Structure, Conductivity, and Ionic Motion in Na_1+xSi_xZr₂P_3−xO₁₂: A Simulation Study, J. Phys. Chem. B, 2002, 106, 7081–7089.
62. S. Roy and P. Padma Kumar, Framework flexibility of sodium zirconium phosphate: role of disorder, and polyhedral distortions from Monte Carlo investigation, J. Mater. Sci., 2012, 47, 4946–4954.
63. S. Roy and P. P. Kumar, Influence of Cationic ordering on ion transport in NASICONs: molecular dynamics study, Solid State Ionics, 2013, 253, 217–222.
64. Z. Deng, G. Sai Gautam, S. K. Kolli, J.-N. Chotard, A. K. Cheetham, C. Masquelier and P. Canepa, Phase Behavior in Rhombohedral NaSiCON Electrolytes and Electrodes, Chem. Mater., 2020, 32, 7908–7920.
65. Z. Zou, N. Ma, A. Wang, Y. Ran, T. Song, Y. Jiao, J. Liu, H. Zhou, W. Shi, B. He, D. Wang, Y. Li, M. Avdeev and S. Shi, Relationships Between Na⁺ Distribution, Concerted Migration, and Diffusion Properties in Rhombohedral NASICON, Adv. Energy Mater., 2020, 10, 2001486.
66. D. T. Qui, J. J. Capponi, J. C. Joubert and R. D. Shannon, Crystal structure and ionic conductivity in Na₄Zr₂Si₃O₁₂, J. Solid State Chem., 1981, 39, 219–229.
67. C. Delmas, R. Olazcuaga, G. Le Flem, P. Hagenmuller, F. Cherkaoui and R. Brochu, Crystal chemistry of the Na_1+xZr_2−xL_x(PO₄)₃ (L = Cr, In, Yb) solid solutions, Mater. Res. Bull., 1981, 16, 285–290.
68. X. He, Y. Zhu and Y. Mo, Origin of fast ion diffusion in super-ionic conductors, Nat. Commun., 2017, 8, 15893.
69. Y. Noda, K. Nakano, H. Takeda, M. Kotobuki, L. Lu and M. Nakayama, Computational and Experimental Investigation of the Electrochemical Stability and Li-Ion Conduction Mechanism of LiZr₂(PO₄)₃, Chem. Mater., 2017, 29, 8983–8991.
70. Y. Noda, K. Nakano, M. Otake, R. Kobayashi, M. Kotobuki, L. Lu and M. Nakayama, Research Update: Ca doping effect on the Li-ion conductivity in NASICON-type solid electrolyte LiZr₂(PO₄)₃: a first-principles molecular dynamics study, APL Mater., 2018, 6, 60702.
71. S. Shi, P. Lu, Z. Liu, Y. Qi, L. G. Hector JR, H. Li and S. J. Harris, Direct Calculation of Li-Ion Transport in the Solid Electrolyte Interphase, J. Am. Chem. Soc., 2012, 134, 15476–15487.
72. S. Roy and P. Padma Kumar, Influence of Si/P ordering on Na⁺ transport in NASICONs, Phys. Chem. Chem. Phys., 2013, 15, 4965–4969.
73. Z. Zhang, Z. Zou, K. Kaup, R. Xiao, S. Shi, M. Avdeev, Y.-S. Hu, D. Wang, B. He, H. Li, X. Huang, L. F. Nazar and L. Chen, Correlated Migration Invokes Higher Na⁺-Ion Conductivity in NaSICON-Type Solid Electrolytes, Adv. Energy Mater., 2019, 9, 1902373.
74. Z. Zou, N. Ma, A. Wang, Y. Ran, T. Song, B. He, A. Ye, P. Mi, L. Zhang, H. Zhou, Y. Jiao, J. Liu, D. Wang, Y. Li, M. Avdeev and S. Shi, Identifying Migration Channels and Bottlenecks in Monoclinic NASICON-Type Solid Electrolytes with Hierarchical Ion-Transport Algorithms, Adv. Funct. Mater., 2021, 31, 2107747.
75. H. Kohler, H. Schulz and O. Melnikov, Composition and conduction mechanism of the NASICON structure X-ray diffraction study on two crystals at different temperatures, Mater. Res. Bull., 1983, 18, 1143–1152.
76. H. Kohler and H. Schulz, NASICON solid electrolytes part I: the Na⁺-diffusion path and its relation to the structure, Mater. Res. Bull., 1985, 20, 1461–1471.
77. H. Kohler and H. Schulz, NASICON solid electrolytes part II – X-ray diffraction experiments on sodium-zirconium-phosphate single crystals at 295 K and at 993 K, Mater. Res. Bull., 1986, 21, 23–31.
78. E. R. Losilla, M. A. G. Aranda, S. Bruque, M. A. París, J. Sanz and A. R. West, Understanding Na Mobility in NASICON Materials: a Rietveld, 23Na and 31P MAS NMR, and Impedance Study, Chem. Mater., 1998, 10, 665–673.
79. H. Park, K. Jung, M. Nezafati, C.-S. Kim and B. Kang, Sodium Ion Diffusion in Nasicon (Na₃Zr₂Si₂PO₁₂) Solid Electrolytes: Effects of Excess Sodium, ACS Appl. Mater. Interfaces, 2016, 8, 27814–27824.
80. S. Naqash, D. Sebold, F. Tietz and O. Guillon, Microstructure–conductivity relationship of Na₃Zr₂(SiO₄)₂(PO₄) ceramics, J. Am. Ceram. Soc., 2019, 102, 1057–1070.
81. Q. Ma, M. Guin, S. Naqash, C.-L. Tsai, F. Tietz and O. Guillon, Scandium-substituted Na₃Zr₂(SiO₄)₂(PO₄) prepared by a solution-assisted solid-state reaction method as sodium-ion conductors, Chem. Mater., 2016, 28, 4821–4828.
82. S. Lunghammer, D. Prutsch, S. Breuer, D. Rettenwander, I. Hanzu, Q. Ma, F. Tietz and H. M. Wilkening, Fast Na ion transport triggered by rapid ion exchange on local length scales, Sci. Rep., 2018, 8, 1–8.
83. M. Kaus, M. Guin, M. Yavuz, M. Knapp, F. Tietz, O. Guillon, H. Ehrenberg and S. Indris, Fast Na⁺ ion conduction in NASICON-type Na_3.4Sc₂(SiO₄)_0.4(PO₄)_2.6 observed by ²³Na NMR relaxometry, J. Phys. Chem. C, 2017, 121, 1449–1454.
84. M. Guin, F. Tietz and O. Guillon, New promising NASICON material as solid electrolyte for sodium-ion batteries: correlation between composition, crystal structure and ionic conductivity of Na_3+xSc₂Si_xP_3−xO₁₂, Solid State Ionics, 2016, 293, 18–26.
85. M. Guin, E. Dashjav, C. M. Kumar, F. Tietz and O. Guillon, Investigation of crystal structure and ionic transport in a scandium-based NASICON material by neutron powder diffraction, Solid State Sci., 2017, 67, 30–36.
86. S. K. Pal, R. Saha, G. V. Kumar and S. Omar, Designing high ionic conducting NASICON-type Na₃Zr₂Si₂PO₁₂ solid-electrolytes for Na-ion batteries, J. Phys. Chem. C, 2020, 124, 9161–9169.
87. J. M. Winand, A. Rulmont and P. Tarte, Ionic conductivity of the Na_1+xM_x^IIIZr_2−x(PO₄)₃ systems (M = Al, Ga, Cr, Fe, Sc, In, Y, Yb), J. Mater. Sci., 1990, 25, 4008–4013.
88. P. J. Squattrito, P. R. Rudolf, P. G. Hinson, A. Clearfield, K. Volin and J. D. Jorgensen, Sodium and oxygen disorder in a scandium-substituted nasicon: a time of flight neutron powder diffraction study of Na_2.5Zr_1.8Sc_0.2Si_1.3P_1.7O₁₂, Solid State Ionics, 1988, 31, 31–40.
89. G. Collin, R. Comes, J. P. Boilot and P. Colomban, Disorder of tetrahedra in Nasicon-type structure—I.: Na₃Sc₂(PO₄)₃: structures and ion-ion correlations, J. Phys. Chem. Solids, 1986, 47, 843–854.
90. M. A. Subramanian, P. R. Rudolf and A. Clearfield, The preparation, structure, and conductivity of scandium-substituted NASICONs, J. Solid State Chem., 1985, 60, 172–181.
91. S. Susman, C. J. Delbecq, T. O. Brun and E. Prince, Fast ion transport in the NASICON analog Na₃Sc₂(PO₄)₃: structure and conductivity, Solid State Ionics, 1983, 9, 839–844.
92. L. Ladenstein, S. Lunghammer, E. Y. Wang, L. J. Miara, H. M. R. Wilkening, G. J. Redhammer and D. Rettenwander, On the dependence of ionic transport on crystal orientation in NASICON-type solid electrolytes, JPhys Energy, 2020, 2, 35003.
93. F. Sun, Y. Xiang, Q. Sun, G. Zhong, M. N. Banis, Y. Liu, R. Li, R. Fu, M. Zheng and T.-K. Sham, Origin of high ionic conductivity of Sc-doped sodium-rich NASICON solid-state electrolytes, Adv. Funct. Mater., 2021, 31, 2102129.
94. T. Zinkevich, A. Fiedler, M. Guin, F. Tietz, O. Guillon, H. Ehrenberg and S. Indris, Na⁺ ion mobility in Na_3+xSc₂(SiO₄)_x(PO₄)_3−x (0.1 < x < 0.8) observed by 23Na NMR spectroscopy, Solid State Ionics, 2020, 348, 115277.
95. Y. Gao, A. M. Nolan, P. Du, Y. Wu, C. Yang, Q. Chen, Y. Mo and S.-H. Bo, Classical and Emerging Characterization Techniques for Investigation of Ion Transport Mechanisms in Crystalline Fast Ionic Conductors, Chem. Rev., 2020, 120, 5954–6008.
96. A. van der Ven, Z. Deng, S. Banerjee and S. P. Ong, Rechargeable Alkali-Ion Battery Materials: Theory and Computation, Chem. Rev., 2020, 120, 6977–7019.
97. X. He, Y. Zhu, A. Epstein and Y. Mo, Statistical variances of diffusional properties from ab initio molecular dynamics simulations, npj Comput. Mater., 2018, 4, 18.
98. R. Kobayashi, nap: a molecular dynamics package with parameter-optimization programs for classical and machine-learning potentials, J. Open Source Softw., 2021, 6, 2768.
99. R. Kobayashi, Y. Miyaji, K. Nakano and M. Nakayama, High-throughput production of force-fields for solid-state electrolyte materials, APL Mater., 2020, 8, 81111.
100. R. Jalem, K. Kanamori, I. Takeuchi, M. Nakayama, H. Yamasaki and T. Saito, Bayesian-Driven First-Principles Calculations for Accelerating Exploration of Fast Ion Conductors for Rechargeable Battery Application, Sci. Rep., 2018, 8, 5845.
101. T. Atsumi, K. Sato, Y. Yamaguchi, M. Hamaie, R. Yasuda, N. Tanibata, H. Takeda, M. Nakayama, M. Karasuyama and I. Takeuchi, Chemical Composition Data-Driven Machine-Learning Prediction for Phase Stability and Materials Properties of Inorganic Crystalline Solids, Phys. Status Solidi B, 2022, 259, 2100525.
102. S. Sekine, T. Hosaka, H. Maejima, R. Tatara, M. Nakayama and S. Komaba, Na[Mn_0.36Ni_0.44Ti_0.15Fe_0.05]O₂ predicted via machine learning for high energy Na-ion batteries, J. Mater. Chem., 2024, 12, 31103–31107.
103. Q. Zhao, L. Zhang, B. He, A. Ye, M. Avdeev, L. Chen and S. Shi, Identifying descriptors for Li⁺ conduction in cubic Li-argyrodites via hierarchically encoding crystal structure and inferring causality, Energy Storage Mater., 2021, 40, 386–393.
104. Q. Zhao, M. Avdeev, L. Chen and S. Shi, Machine learning prediction of activation energy in cubic Li-argyrodites with hierarchically encoding crystal structure-based (HECS) descriptors, Sci. Bull., 2021, 66, 1401–1408.
105. S. Takamoto, C. Shinagawa, D. Motoki, K. Nakago, W. Li, I. Kurata, T. Watanabe, Y. Yayama, H. Iriguchi, Y. Asano, T. Onodera, T. Ishii, T. Kudo, H. Ono, R. Sawada, R. Ishitani, M. Ong, T. Yamaguchi, T. Kataoka, A. Hayashi, N. Charoenphakdee and T. Ibuka, Towards universal neural network potential for material discovery applicable to arbitrary combination of 45 elements, Nat. Commun., 2022, 13, 2991.
106. R. Iwasaki, N. Tanibata, H. Takeda and M. Nakayama, Universal-neural-network-potential molecular dynamics for lithium metal and garnet-type solid electrolyte interface, Commun. Mater., 2024, 5, 148.
107. J. Qi, T. W. Ko, B. C. Wood, T. A. Pham and S. P. Ong, Robust training of machine learning interatomic potentials with dimensionality reduction and stratified sampling, npj Comput. Mater., 2024, 10, 43.
108. C. Chen and S. P. Ong, A universal graph deep learning interatomic potential for the periodic table, Nat. Comput. Sci., 2022, 2, 718–728.
109. K. Makino, N. Tanibata, H. Takeda and M. Nakayama, Computational studies on Mg ion conductivity in Mg_2xHf_1−xNb(PO₄)₃ using neural network potential, J. Solid State Electrochem., 2024, 28, 4471–4481.
110. D. Kwon and D. Kim, Machine learning interatomic potentials in engineering perspective for developing cathode materials, J. Mater. Chem., 2024, 12, 23837–23847.
111. A. Hjorth Larsen, J. Jørgen Mortensen, J. Blomqvist, I. E. Castelli, R. Christensen, M. Dułak, J. Friis, M. N. Groves, B. Hammer, C. Hargus, E. D. Hermes, P. C. Jennings, P. Bjerre Jensen, J. Kermode, J. R. Kitchin, E. Leonhard Kolsbjerg, J. Kubal, K. Kaasbjerg, S. Lysgaard, J. Bergmann Maronsson, T. Maxson, T. Olsen, L. Pastewka, A. Peterson, C. Rostgaard, J. Schiøtz, O. Schütt, M. Strange, K. S. Thygesen, T. Vegge, L. Vilhelmsen, M. Walter, Z. Zeng and K. W. Jacobsen, The atomic simulation environment—a Python library for working with atoms, J. Phys.:Condens. Matter, 2017, 29, 273002.
112. H. Chen, L. L. Wong and S. Adams, SoftBV–a software tool for screening the materials genome of inorganic fast ion conductors, Acta Crystallogr., Sect. B, 2019, 75, 18–33.
113. H. Chen and S. Adams, Bond softness sensitive bond-valence parameters for crystal structure plausibility tests, IUCrJ, 2017, 4, 614–625.
114. F. H. Stillinger and T. A. Weber, Computer simulation of local order in condensed phases of silicon, Phys. Rev. B:Condens. Matter Mater. Phys., 1985, 31, 5262.
115. K. Nakano, N. Tanibata, H. Takeda, R. Kobayashi, M. Nakayama and N. Watanabe, Molecular dynamics simulation of Li-ion conduction at grain boundaries in NASICON-Type LiZr₂(PO₄)₃ solid electrolytes, J. Phys. Chem. C, 2021, 125, 23604–23612.
116. Z. Yang, R. E. Ward, N. Tanibata, H. Takeda, M. Nakayama and R. Kobayashi, Exploring the diffusion mechanism of Li ions in different modulated arrangements of La_(1−x)/3Li_xNbO₃ with fitted force fields obtained via a metaheuristic algorithm, Solid State Ionics, 2021, 366, 115662.
117. S. Miyakawa, S. Matsuda, N. Tanibata, H. Takeda, M. Nakayama, T. Saito and S. Fukuchi, Computational studies on defect chemistry and Li-ion conductivity of spinel-type LiAl₅O₈ as coating material for Li-metal electrode, Sci. Rep., 2022, 12, 16672.
118. S. Aizu, S. Takimoto, N. Tanibata, H. Takeda, M. Nakayama and R. Kobayashi, Screening chloride Li-ion conductors using high-throughput force-field molecular dynamics, J. Am. Ceram. Soc., 2023, 106, 3035–3044.
119. X.-S. Yang and S. Deb, in 2009 World congress on nature & biologically inspired computing (NaBIC), IEEE, pp. , pp. 210–214.
120. X.-S. Yang and S. Deb, Engineering optimisation by cuckoo search, Int. J. Math. Model. Numer. Optim., 2010, 1, 330–343.
121. X.-S. Yang and S. Deb, Cuckoo search: recent advances and applications, Neural Comput. Appl., 2014, 24, 169–174.
122. J. Rodrigues and A. Zellner, Weighted balanced loss function and estimation of the mean time to failure, Commun. Stat. Theor. Methods, 1994, 23, 3609–3616.
123. P. Hohenberg and W. Kohn, Inhomogeneous electron gas, Phys. Rev., 1964, 136, B864.
124. W. Kohn and L. J. Sham, Self-consistent equations including exchange and correlation effects, Phys. Rev., 1965, 140, A1133.
125. G. Kresse and J. Hafner, Ab initio molecular dynamics for liquid metals, Phys. Rev. B:Condens. Matter Mater. Phys., 1993, 47, 558.
126. G. Kresse and J. Hafner, Norm-conserving and ultrasoft pseudopotentials for first-row and transition elements, J. Phys.:Condens. Matter, 1994, 6, 8245.
127. G. Kresse and J. Furthmüller, Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set, Phys. Rev. B:Condens. Matter Mater. Phys., 1996, 54, 11169.
128. G. Kresse and J. Furthmüller, Efficiency of ab initio total energy calculations for metals and semiconductors using a plane-wave basis set, Comput. Mater. Sci., 1996, 6, 15–50.
129. M. Parrinello and A. Rahman, Crystal Structure and Pair Potentials: A Molecular-Dynamics Study, Phys. Rev. Lett., 1980, 45, 1196–1199.
130. M. Parrinello and A. Rahman, Polymorphic transitions in single crystals: a new molecular dynamics method, J. Appl. Phys., 1981, 52, 7182–7190.
131. J. P. Perdew, K. Burke and M. Ernzerhof, Generalized gradient approximation made simple, Phys. Rev. Lett., 1996, 77, 3865.
132. P. E. Blöchl, Projector augmented-wave method, Phys. Rev. B:Condens. Matter Mater. Phys., 1994, 50, 17953–17979.
133. G. Kresse and D. Joubert, From ultrasoft pseudopotentials to the projector augmented-wave method, Phys. Rev. B:Condens. Matter Mater. Phys., 1999, 59, 1758.
134. H. J. Monkhorst and J. D. Pack, Special points for Brillouin-zone integrations, Phys. Rev. B, 1976, 13, 5188.
135. J. Schuett, A. S. Kuhn and S. Neitzel-Grieshammer, Predicting the Na⁺ ion transport properties of NaSICON materials using density functional theory and Kinetic Monte Carlo, J. Mater. Chem., 2023, 11, 9160–9177.
136. J. Schuett, F. Pescher and S. Neitzel-Grieshammer, The origin of high Na⁺ ion conductivity in Na_1+xZr₂Si_xP_3−xO₁₂ NASICON materials, Phys. Chem. Chem. Phys., 2022, 24, 22154–22167.
137. Y. Deng, C. Eames, L. H. B. Nguyen, O. Pecher, K. J. Griffith, M. Courty, B. Fleutot, J.-N. Chotard, C. P. Grey and M. S. Islam, Crystal structures, local atomic environments, and ion diffusion mechanisms of scandium-substituted sodium superionic conductor (NASICON) solid electrolytes, Chem. Mater., 2018, 30, 2618–2630.
138. A. Pedone, G. Malavasi, M. C. Menziani, A. N. Cormack and U. Segre, A new self-consistent empirical interatomic potential model for oxides, silicates, and silica-based glasses, J. Phys. Chem. B, 2006, 110, 11780–11795.
139. A. L. Usler, D. Kemp, A. Bonkowski and R. A. de Souza, A general expression for the statistical error in a diffusion coefficient obtained from a solid-state molecular-dynamics simulation, J. Comput. Chem., 2023, 44, 1347–1359.
140. A. Bito, in IEEE Power Engineering Society General Meeting, 2005, vol. 2, pp. 1232–1235.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d5ra01549a

---