Ring Mechanism of Fast Na Ion Transport in Na2LiFeTeO6: Insight from Molecular Dynamics Simulation

Kartik Sau 2, 3, ∗ and Tamio Ikeshoji Research Institute of Electrochemical Energy (RIECEN), National Institute of Advanced Industrial Science and Technology (AIST), 1-8-31 Midorigaoka, Ikeda, Osaka 563-8577, Japan Mathematics for Advanced Materials Open Innovation Laboratory (MathAM-OIL), National Institute of Advanced Industrial Science and Technology (AIST), c/o Advanced Institute for Materials Research (AIMR), Tohoku University, Sendai 980-8577, Japan Advanced Institute for Materials Research (AIMR), Tohoku University, Sendai 980-8577, Japan


I. INTRODUCTION
The improvement of high-energy-density batteries is a major focus in current research [1] owing to consumer demand for sophisticated electronic gadgets and the everincreasing demand for electric vehicles. Conventional Liion batteries cannot satisfy present demand, whereas allsolid-state batteries(ASSB) could be suitable because of favorable factors such as high volumetric and gravimetric energy density, safety, and long life cycle. Fast ion conducting (FIC) solids form an integral part of ASSB and are a leading topic of current research [2][3][4]. Significant efforts have been devoted over the last three decades to the advancement of Li + ion conducting materials [5][6][7][8][9][10][11][12][13][14]. Recently, sodium ion conductors have been reported to be promising for energy storage applications because the Na ion has a very similar intercalation chemistry to the Li ion and, in some cases, the Na ionic radius is more suitable for stabilizing fast ion conducting frameworks because of its larger size compared with the Li ionic radius. Thus, the development of sodium ion conducting ASSB working at ambient temperatures is gaining renewed attention.
Recently, a unique class of heterostructures known as honeycomb layered oxides has gained attention as highenergy-density electrodes and electrolyte materials for their rich crystal chemistry that provides high cationic conduction and high voltages using Na and K ions [15][16][17][18][19][20][21][22][23][24][25][26][27]. The general formula of this series of materials is A 2 M 2 DO 6 , A 3 M 2 DO 6 or A 4 M DO 6 , where M can be divalent or trivalent transition metal atoms such as Cr, Mn, Fe, Co, Ni, Cu, or some combination thereof; D rep- * kartik.sau@gmail.com resents pentavalent or hexavalent metal atoms such as Sb, Te, or Bi; and A can be alkali metal atoms such as Li, Na, or K [15,28]. The structure in general consists of an array of parallel transition metal slabs with a distinct honeycomb arrangement of multiple M atoms surrounding D atoms in a layered framework of interposed A alkali atoms. Specifically, the tellurates in honeycomb layered oxides have been reported to produce the highest voltages (over 4V) to date [16-18, 20-23, 26]. Thus, a thorough study on the cationic diffusion mechanism in honeycomb layered oxides based on tellurate (Na 2 LiFeTeO 6 ) presents an ideal platform for maximizing the cationic conduction.
Atomistic details are valuable for better understanding and development of FIC solids [29]. Molecular dynamics (MD) simulation has emerged as one of the first simulation methods for its ability to elucidate the underlying mechanisms at atomic level [30], as a complement to experimental probes, which are often difficult to implements. Silver iodide [31,32], β-alumina [33], and Na 1+x Zr 2 Si x P 3−x O 12 (NASICON) (0 ≤ x ≤ 3) [34] are classic examples that show how MD studies predicted the experimental results. The technique has been successfully applied to a variety of oxide ion conductors as well [35,36]. Recently, we developed an inter-atomic potential model for a family of fast ion conductors-A 2 M 2 TeO 6 (A = Li, Na, or K ; M = Ni, Zn, Co, or Mg) [37,38]-to perform MD studies that represent the structural and dynamical properties reported in the experiments and present several interesting behaviors such as energy-driven or entropy-driven diffusion depending on the concentration of the mobile ions. Our previous work on A 2 Ni 2 TeO 6 (A = Li, Na, or K) [38] showed how cationic distribution and the energy landscape vary with increasing inter-layer distance. The study also showed that ionic conductivity increases exponentially with increasing inter-layer distance. Moreover, the simulated cationic path in Na 2 Ni 2 TeO 6 has been experimentally confirmed [39,40]. Accurately depicting microscopic physicochemical properties of materials that are beyond experimental reach requires the use of reliable potential parameters that can effectively reproduce the structural and transport properties determined by experiment. In this context, first-principles MD (FPMD) simulation, which reproduce forces on atoms precisely from electronic structure, plays a major role. However, FPMD is computationally extremely expensive and cannot deal with larger system sizes. Therefore, achieving accurate diffusion behavior is difficult because of low statistical accuracy. Thus, we rely on force-field-based MD simulations.
The present MD simulations are motivated by the recent experimental studies by Nalbandyan et al. [41], who synthesized a fast Na + -conducting solid (0.04 S/cm at 573 K comparable with the β-alumina [42][43][44][45]) of the formula Na 2 LiFeTeO 6 . The material is isostructural with heterovalence analogs of A 2 M 2 TeO 6 (A = Na, K ; M = Ni, Zn, Co, and Mg), consisting of edge-sharing LiO 6 , FeO 6 and TeO 6 octahedral brucite-like layers spanning the ab-plane, in which Te alternates with Li or Fe. The brucite-like octahedral layers are not identical along the c-axis, unlike the structure A 2 Ni 2 TeO 6 (A = Na, K), where the Na + ions occupy the inter-layer. The present study provides full inter-atomic potentials and fresh insights on ion transport.
To illustrate the local octahedral ordering, a crystal structural illustration of Na 2 LiFeTeO 6 is furnished in figure 1. As shown, the Na + ions are interspersed between metal slabs consisting of Li, Fe, and Te in an octahedral coordination environment with oxygen atoms. Despite the immense potential envisioned in the experimental study, details theoretical explorations of honeycomb layered materials, such as surrounding of the various mesoscopic mechanisms engendering their remarkable electrochemical performance still remain unexplored. Specifically, local octahedral inter-layer influencing on ionic conductivity are unclear, as the inter-layer distance is determined by existing elements in the framework octahedral layer and by mobile ions within the framework inter-layer.

A. Framework Polyhedral Structure
The Na 2 LiFeTeO 6 structure consists of several polyhedral inter-layers with sandwiching Na + ions, as shown in Figure 1 reported by Nalbandyan et al. [41]. The top and the bottom of the adjacent framework inter-layers are not identical, unlike the other materials, A 2 Ni 2 TeO 6 (A = Na, and K). The LiO 6 , FeO 6 , and TeO 6 octahedra are arranged in a honeycomb fashion inside the framework layer. The octahedral inter-layers are stabilized along the c-axis by Van der Waals forces and by interactions me- . Blue lines denote the unit cell and the inter-layer distance is marked by d. (b) Several crystallographically distinct sites with face-sharing framework octahedra in which Na ions reside, denoted as Na1 (in orange), Na2 (in yellow), and Na3 (in purple). (c) Na sites Na4 (in magenta), Na5 (in green), and Na6 (in cyan), with edge-sharing framework octahedra. The local polyhedral environments are also shown, for clarity. The a-and c-axes are also shown.
diated via the Na atoms occupying the inter-layers. To understand octahedral ordering, an octahedral density pattern of Na 2 LiFeTeO 6 along the y-direction from MD simulation is shown in Figure 2, in which the FeO 6 octahedra are located on the top and the bottom along the c-axis. The Na + ion density pattern in the conduction plane is significantly influenced by the local octahedral sequences, most importantly where the inter-layer distance is small (∼ 3.5Å). The octahedral density pattern from MD simulation follows same ordering as reported in X-ray structure. The inter-layer octahedral density are found to be a slightly off particularly at high temperatures because of thermal oscillation that leads to polyhedral layer-sliding. Furthermore, the average cell parameters of Na 2 LiFeTeO 6 calculated from NPT-MD simulations at 300 K are comparable to the experimental reported cell parameters, with about 1.5% error, as shown in Table I, validating the potential model. To obtain a more detailed structural comparison, we also calculated radial distribution functions, g(r), between the frame- Atomic distribution of octahedral center atoms along the y-axis. The octahedral center atoms-Li, Te, and Fe-are marked in light green, blue, and brown, respectively, at 300 K for Na 2 LiFeTeO 6 from NPT-MD force-field-based simulations. The corresponding crystallographic octahedral layers [41] are also shown in ball-and-stick model using the same color coding for comparison. % 800 K for Na 2 LiFeTeO 6 . All calculated peaks from MD simulations are consistent with their corresponding g(r), calculated for the XRD structure and FPMD, although thermal broadening is observed. Whereas, a sharper peak from force-field-based MD simulations indicates stiffer inter-atomic potential parameters. It could be possibly because of the structural disordered (different B-O bond distance; B = Li, Fe, and Te), whereas we employed a single type of pair of interaction for simplicity of the potential model. The Li, Fe, and Te atoms have exactly six O atoms within their respective bond distances, as shown in Figure S1 in the Supplementary Information. The bond lengths are also close to the experimental reported bond distances, affirming that the present set of inter-atomic parameters reasonably represents the local structure as well as the octahedral sequences. Herein, we want to emphasize that most of the cases the Li-atoms are consider to be mobile; in this case they form immobile framework. In contrast, the loose coordination to the framework polyhedral inter-layer allocates several cationic sub-lattices in the conduction layer for facile cationic conduction, as identified in the reported experimental structure. The facile Na + ion diffusion is reflected in the plot of the mean square displacements (MSD) against time (figure 4(a)). The closely packed framework octahedral layers parallel to the ab-plane restrict the Na + ion diffusion along the c-axis direction of the cell, as reflected in the inset of Figure 4(b) (less than 0.04Å 2 [46]). Therefore, cationic mobility appears to be restricted within the sublattices oriented parallel to the ab-plane, affirming the diffusion to be highly anisotropic (or layer diffusion), as also observed in the previous studies [20,47].
The Na + ion conductivity, σ, is calculated from the diffusion coefficient, D (the slope of the MSD near the diffusive region), according to Eq. 4. Calculation of the conductivity is also possible using the Green-Kubo relation instead of the Nernst-Einstein eq. 3. In principle, both approaches are equivalent. However, the correlation function involved in the Green-Kubo formula shows Inset shows the MSD of Na + ions along the c-axis direction (MSD-z) at 600 K. These results demonstrate that the framework atoms are immobile and the Na + ions are restricted along the ab-plane, rendering the diffusion to be highly anisotropic.
a slow decay or what is referred to as long-time tail [48] that includes large fluctuations in the calculated electrical conductivity at long times [49][50][51]. The plot of the logarithm of σT versus the inverse temperature (1000/T ) is displayed in Figure 5. The simulated result displays a straight line indicating Arrhenius behavior, whereas the experimental result shows a deviation from linearity at high temperature, possibly because of a phase transition. The activation energy can be calculated from the slope of the straight line according to Eq. 3. The activation energy of Na 2 LiFeTeO 6 is found to be 0.54 ± 0.02 eV from the MD simulation, whereas 0.48 eV is reported from the experimental study [41]. Therefore, we conclude that the present MD simulation model successfully represents the structural (framework octahedral orientation and radial distribution functions) and dynamical properties of Na 2 LiFeTeO 6 with reasonable experimental agreement.

C. Atomistic Mechanisms of Na + Ion Transport
Herein, our focus is on the atomistic origin of high ionic conductivity in Na 2 LiFeTeO 6 . We calculated the threedimensional isosurface population density along with the polyhedral framework, as shown in Figures 6 (a) and (b) for understanding of Na + ion transport path. The Na + ion density lies inside the octahedral framework layer, which is also reflected in the MSD plot of the Na + ions along the z-direction. For further clarity of the Na + ions' diffusion paths, the two-dimensional population density (the number of ions per unit area) pattern and the average potential energy of the Na + ions in the conduc-  The octahedral framework structure is also shown. (c) Average population density profiles of the Na + ion for Na 2 LiFeTeO 6 in the ab-plane derived from NVE-MD simulations at 600 K (mapped onto 2×2 unit cells). The locations of the crystallographically distinct Na sites, namely Na1, Na2, Na3, Na4, Na5, and Na6, are shown in cyan. Three high-density sites identified in MD are marked by S1, S2, and S3. The locations of the crystallographic octahedral centers of Li (orange), Fe (green), and Te (red) are also marked by filled circles. The population density profile indicates two preferred migration paths of the Na + ions: (1) intra-ring, consisting of two (ring in orange) or three (ring in sky blue) high-density areas and (2) inter-ring, moving from the orange ring to the sky blue ring. The high-density sites S2 and S3 are located inside the purple triangles formed by connecting the octahedral centers, whereas the S1 site is located on an edge of a triangle. (d) Corresponding average potential energy (eV) profile of Na + ions in the ab-plane. The Na + ion avoids occupying potentially favorable sites by reflecting a low population at potential minima. The locations of the high-density populated sites (S1, S2, and S3) are also marked by circles, as shown in the legends. The unit cell is marked by dotted lines.
The average potential energy of the individual cations is calculated using Eq. 6 (detailed in the METHODOL-OGY section IV) and projected onto a 2D grid in the ab-plane. Both the projected population density and the average potential energy are also replicated in 2×2 unit cells in the conduction plane for continuity, as shown in Figures 6(c) and (d). The high density areas and their local environment-the octahedral centers and the Na sites-are indicated on top of the density and energy patterns. There are five high-density maxima in the population density profile, whereas there are four available Na + ions per unit cell per inter-layer. Therefore, the Na + ions can occupy the high-density maxima and can diffuse as well through the available vacant sites; as a result, the occupancy at each high-density area has to be less than unity. The occupancy at each maximaum density site was calculated using Bader charge analysis code, which was principlly developed for electron density analysis [52]. Herein, the iso-surface density of Na + ion was used instead of charge density. The occupancy at S1 and S2 are close to unity, whereas close to half occupancy is identified at S3. For cationic channel connectivity, two types of ring connectivity are observed, as marked by circles: ring 1 (sky blue circle) consists of three high-density maxima, whereas ring 2 (orange circle) is formed by two maxima. The sites are marked S1, S2, and S3, as shown in Figure 6(c). Among the three high-density sites, S1 and S3 find good agreement with the reported crystallographic sites, whereas site S2 deviates from any reported crystallographic site. Furthermore, the sites S2 and S3 are located inside the circles formed by the framework octahedral centers, whereas site S1 is identified on the side of a triangle. None of the Na sites are found on top of the octahedral center because of strong repulsion between the top and bottom octahedra and the Na + ions.
The high-density population maxima are indicated on top of the potential energy profile in Figure 6(d), which reveals a different behavior than the population profile; the potential energy minima do not match with the population maxima. Rather, the high densities are found around the potential energy minima, possibly because of the entropic effect. For instance, two potential energy minima per layer per unit cell are identified, whereas four Na + ions are available. Thus, the sodium ions cannot occupy the energetically favorable sharp potential energy regions (-2.4 eV), but rather occupy the sites which are greater in number and therefore higher in entropy. The difference between the effects of energy and entropy can be linked through Eq. 8, in which entropy plays a significant role. Specifically, particles usually avoid energetically favorable sharp potential energy minima with low entropy, whereas they prefer a relatively higher energy basin which has high entropy. This scenario can be understood by implementing a simulation which places fewer Na + ions in a conducting layer (detailed in the METHODOLOGY section IV). The Na + ion density pattern and the potential energy profile are shown in Figure  S2 in the supplementary information for the under-loaded TABLE II. Free energy barrier of Na + ions hopping between two population density maxima, as marked on Figure 6(c), calculated with Eq. 9 at 600 K.
Between the Sites Site Distance (Å) Barrier Height (eV) Ring 1 S2-S3 2.7 0.23 S3-S3 2.2 0.14 Ring 2 S1-S1 3.0 0.25 Inter Ring S2-S1 3.2 0.34 conduction plane. In the potential energy profile, the locations of the potential energy minima are almost unchanged, whereas the population density pattern changes significantly. Only one type of ring exists, consisting of four high-density maxima. The highest population densities are identified exactly at the potential minima. In this case, the smaller number of Na + ions contributes to a lower entropy, resulting in an energy-driven population distribution. The entropy contribution is even lower if the number of Na + ions in the conduction layer is reduced even further, as shown in Figure S3, where all population density maxima and potential minima coincide exactly in the conduction plane, implying a minimal entropic contribution. For a quantitative estimation of the intra-and interring free energy barrier heights that determine the Na + ion diffusion mechanism, the one-dimensional free energy distributions among the high-density population maxima are shown in Figure 7. The free energies were calculated by counting the population density distributions inside a rectangle connecting two high-density areas (with a width of 1.0Å) as indicated in Figure 6 and projected along a straight line connecting two sites using Eq. 9. Figure 7 captures the free energy minima and maxima for inter and intra-ring cases, distinctly, as reflected in Figure 6. The pairs of site distances and corresponding free energy barriers of Na + ion hopping are listed in Table II. The inter-ring free energy barriers are higher than the intra-ring free energy barriers, which is also reflected in Figure 6. Thus, the intra-ring mechanism is more favorable than the inter-ring mechanism. However, both mechanisms are responsible for the long-range Na + ion diffusion. A direct ring mechanism is also revealed by randomly tracking a Na + ion over simulation time, as displayed in Figure 8, where the high-density Na sites (S1, S2, and S3) are also marked. The Na + ions diffuse through both the inter-ring and intra-ring pathways, as evident in the free energy barriers and the population density profiles.
Both the population density profile and the free energy profile indicate that the Na + ion is trapped for a while within the ring. Therefore, the estimation of the trapping time of the Na + ions provides a valuable insight into Na + ion diffusion. Thus, the self-van Hove correlation function is calculated (Figure 9(a)). The first peak in this figure indicates the thermal oscillation of Na + ions in their respective lattice sites, whereas the other peaks signify the Na + ions hopping to the neighboring sites. The residence time (the time taken before the second peak arises) of Na + ions at the lattice sites is observed to be roughly 10 ∼ 15 ps. However, the third peak appears after a significantly longer time (∼ 100 ps), because the Na + ions become trapped inside the ring. Furthermore, the trapping time of Na + ions inside the ring is also found in a log-log plot of the Na + ions' MSD, as shown in Figure  9(b), which shows a caging of about 100 ps (roughly until the point at which the log-log MSD follows straight-line behavior) with an MSD value of about 4Å 2 before moving to the diffusive region. Interestingly, this 4Å 2 is close to a value of d 2 /2, where d is the approximate diameter of a ring (∼ 3Å). Therefore, the log-log MSD plot also supports the ring mechanism of the diffusion of the Na + ions. The above behavior is also confirmed in the underloaded Na plane for case 1, as displayed in Figure S3 in the supplementary information. In this case, the Na + ion also shows trapping inside a ring consisting of four high-density maxima, as revealed in Figure S2(a) in the Supplementary Information. The corresponding self-van Hove correlation function also reflects a trapping time of about 50 ∼ 60 ps in Figure S3(a). The same time duration is also identified in the log-log MSD plot in Figure  S3(b).
For further clarity of the Na + ion hopping mechanism, the coordination number of the Na + ion is analyzed. The first peak at 3Å in Figure 10(a) indicates that if one S3 site is occupied, the nearest S3 site has to remain vacant (the nearest S3-S3 distance is 2.25Å). The above scenario can be explained using a schematic diagram, as shown in Trajectory of a randomly chosen Na + ion (purple) on the ab-plane at 600 K. The maximum density MD sites, S1 (red), S2 (green), and S3(cyan), are also marked. The unit cells are indicated by black lines. Figure 10(b). The most likely distribution of Na + ions invokes the occupation of two sites of a ring, leaving one site vacant to avoid strong Na-Na repulsion, as there are six Na sites and four Na + ions per layer per unit cell. Thus, the Na + ions can show a circular motion inside the ring, or they can move out from the ring. When one Na + ion moves out from a ring, another Na + ion should enter from a neighboring ring to maintain local charge  9. (a) The self-van Hove correlation function, Gs(r, t ), as a function of radial distance (r) and time (t ) for Na 2 LiFeTeO 6 at 600 K from NVE-MD simulations, shown on a logarithmic scale for better visualization. Several high-intensity regions (marked as 1st, 2nd, and 3rd) indicate the Na + ion occupancy while hopping from one site to another site occurs. Time-scales of just before grow the 2nd and 3rd peaks are also mentioned. (b) A log-log plot of MSD for identifying the ion trapping in a ring. A dotted red line is added to distinguish the caging and diffusive regions.
balance. If these ions show a circular motion within a ring, the duration should be close to 100 ps, as identified in the log-log plot of the MSD in Figure 9(b).

III. CONCLUSION
A refined set of inter-atomic models represents the structural and transport properties of honeycomb layered oxides in Na 2 LiFeTeO 6 , in excellent agreement with the experimental results. This potential model can be leveraged to garner distinct atomistic insights such as the entropic contribution in cationic distribution and the ring-like cationic diffusion from the population density profile, the free energy barrier, the self-van Hove correlation, and the log-log MSD plot. Therefore, an MD simulation provides a new avenue to test the ion dynamics of various honeycomb layered oxides. Particularly, the entropic contribution and the distinct ring-like feature that controls fast ion transport in solids is interesting in applications beyond honeycomb layered oxides.

A. Interatomic Potential
To carry out an MD simulation, the reported Vashishta-Rahman form of interatomic potential [31] is employed. This method is becoming popular recently, as it describes structural and dynamical properties of a series of honeycomb layered oxides, namely Na 2 M 2 TeO 6 (M = Ni, Zn, Co and Mg), [37] as well as NASICONs.  A Na−Na = 0 eV, A Na−Li = 80000 eV, A Na−Fe = 27000 eV, A Na−Te = 7000 eV, n Na−Na = n Na−Li = n Na−Fe = n Na−Te = 11 a The value used from the reported literature [37,53].
This potential form is softer in nature than the other widely used Born-Mayer (Buckingham) [54] and Lennard-Jones potentials [55] The Vashishta-Rahman potential form is as follows [31]; where q i represents the charge andσ i the ionic radius of the i-th atom. The parameters A ij , P ij and C ij describe the overlap-repulsion energy of the electron clouds, the average charge dipole interactions, and the dispersion constant between the ion pairs i and j, respectively. The Vashishta-Rahman potential has a softer overlap repulsion (1/r n , where n = 11, 9, or 7 for cation-cation, cationanion, and anion-anion pairs, respectively), particularly between the anion-anion pairs. The inter-atomic parameters used in this study are listed in Table III. Some of the parameters have already been reported [37]. The parameters that are unavailable in the literature were determined using empirical fitting to attain the experimentally reported bond lengths (Li-O, Fe-O, and O-O) of Na 2 LiFeTeO 6 [34,37,56,57]. The details of the empirical fitting methods are as follows: The non-identical brucite-like octahedral framework along the c-axis is retained by refining the Li-Na, Fe-Na, and Te-Na pairs and the Na-O parameters were refined to represent the conductivity.

B. Computational Details
The empirically fitted set of parameters listed in Table III were used to carry out a series of MD simulations in the temperature range of 400 K to 600 K with 25 K intervals, and zero atmospheric pressure. The Parrinello-Rahman isobaric-isothermal (NPT) MD method [58], which allows for changes in the simulation box sizes while keeping angles fixed, was used. The temperatures and pressure in the system were controlled using thermostating and barostating techniques, whereby some dynamic variables are coupled with particle velocities and simulation box dimensions. Simulations commenced from the reported X-ray diffraction structures [41] wherein the simulation supercell was constructed by 7 × 4 × 3 unit cells comprising 3696 atoms in an orthorhombic symme-try following P2 1 2 1 2 1 (No. 19). Several partially occupied crystallographic sites were identified in the Na conduction plane. In MD simulations, we placed the Na + ions at the crystallographically reported highest occupied sites (Na2 and Na5) to avoid strong cation-cation repulsion. For clarity, all reported crystallographically distinct sites (Na1, Na2, Na3, Na4, Na5, and Na6) [41] in the layered frameworks of Na 2 LiFeTeO 6 are shown in Figure 1. The inter-layer interaction is usually weaker in nature; thus, inter-layer sliding is expected unless the initial configuration is not set properly. To avoid inter-layer sliding, we performed a constrained MD simulation for 500 ps for each temperature case; the framework atoms (Li, Fe, and Te) were frozen at the crystallographic sites and the other atoms were allowed to move. Finally, we lifted the former constraint and performed the usual MD simulation. A time step of 2 fs was applied to the velocity Verlet algorithm to solve Newton equations of motion. A typical run-time of 6 ns or longer was used with trajectory samples stored at intervals of 200 fs for detailed analyses. To guarantee the thermodynamic convergence properties, a few longer run-time simulations of 100 ns were performed separately. Periodic boundary conditions in all three directions and the Ewald summation technique for the convergence of long-range Coulombic interactions were applied using the LAMMPS software package [59]. A cut-off distance of 13Å was used for both the short-range interactions and the short-range part of the Ewald summation. Micro-canonical MD (NVE-MD) simulations were further performed at 600 K using the final structure obtained from NPT-MD simulations.
Furthermore, to understand the entropy effect in the presence of a larger number of Na + -ions, we performed a few constrained MD simulations with a different number of Na + ions in each layer. For instance, there are four Na + -ions in each layer in each unit cell. Here, we performed two distinct NPT-MD simulations with alternative Na + -ion densities as follows: 1. The first layer contains three Na + ions (normally four), whereas five Na + -ions were placed at Na sites in the next layer. We call this scenario underloaded case 1.
2. The first layer contains two Na + ions (normally four), whereas six Na + -ions were placed at Na sites in the next layer. We call this scenario underloaded case 2.
We constructed two such super-cells additionally, which also maintain charge conservation, and performed NVE-MD simulations followed by NPT-MD simulations at 600 K.
Finally, an FPMD simulation was performed based on Density Functional Theory (DFT) using the Vienna Ab Initio Simulation Package (VASP) [60,61] to compare the structure with that obtained using the force-field-based MD simulation. We performed FPMD simulations in an NPT ensemble, in which temperature and pressure were controlled using Langevin thermostating and barostating techniques. The plane-wave basis sets and projectoraugmented-wave pseudopotentials [62,63] were used under periodic boundary conditions. An energy cutoff of 320 eV was used; integrals over the Brillouin zone were performed only at the Γ-point. The simulated super-cell consisted of 2× 1× 1 unit cells with a total of 88 atoms. The FPMD simulation was performed at 800 K for 1-ps equilibration and 2-ps sampling with a time step of 1 fs.

C. Estimation of Key Quantities
The ion transport is calculated from the mean square displacement (MSD) of Na + following Einstein's relation. The self-diffusion coefficient is calculated from the slope of the MSD for the two-dimensional transport as where N denotes the total number of mobile atoms in the system, r j (t) is the position vector of the j th ion at time t, t is the time difference, and the angular bracket indicates the average over several points in time. The factor four appeared in the denominator of eq. 2b is because of two-dimensional diffusion, whereas it is six for the three-dimensional case. The diffusion coefficient, D, depends on the temperature (T ) according to the Arrhenius equation, where D 0 is the pre-exponential factor, E a represents the activation energy of ion hopping, and k B is the Boltzmann constant. The conductivity, σ, can be linked to the diffusion coefficient, D, according to the Nernst-Einstein equation, where n is the carrier density, q represents the formal charge of the mobile cation, and H r is Heaven ratio, which is assumed to be one for the present system. By substituting these relationships into Eq. 3, one can write ln(σT ) = ln( The potential energy of the individual cations is calculated as where V ij is the interaction potential given in Eq. (1), such that the total potential energy of the system is The relative free energy as seen by the individual cations, ∆F , can be measured relative to another system of cations at the location on the ab-plane exhibiting the maximum cation population density, p max , according to the formula [64][65][66][67], where p represents the population density function (the occupancy) of the cations, p max is the maximum value of the population density function, ∆V is the difference in potential energy of the individual cations (eq. (4)), and ∆S is the difference in entropy (the relative entropy). Finally, we calculated the "self" part of the self-van Hove correlation function, G S (r, t ), to understand the Na + ion hopping mechanism using the formula [68], G S (r, t ) = 1 N N j=1 δ( r − ( r i (t + t)) − r i (t))) . (9)