Modelling Cationic Diffusion in Nickel-Based Honeycomb Layered Tellurates using Vashishta-Rahman Interatomic Potential and Relevant Insights
Although the fascinatingly rich crystal chemistry of honeycomb layered oxides has been accredited as the propelling force behind their remarkable electrochemistry, the atomistic mechanisms surrounding their operations remain unexplored. Thus, herein, we present an extensive molecular dynamics study performed systematically using a reliable set of inter-atomic potential parameters of A2Ni2TeO6 (where A = Li, Na, and K). We demonstrate the effectiveness of the Vashishta-Rahman form of the inter-atomic potential in reproducing various structural and transport properties of this promising class of materials and predict an exponential increase in cationic diffusion with larger inter-layer distances. The simulations demonstrate the correlation between broadened inter-layer (inter-slab) distances associated with the larger ionic radii of K and Na compared to Li and the enhanced cationic conduction exhibited in K 2 Ni 2 TeO 6 and Na 2 Ni 2 TeO 6 relative to Li 2 Ni 2 TeO 6 . Whence, our findings connect lower potential energy barriers, favourable cationic paths and wider bottleneck size along the cationic diffusion channel within frameworks (comprised of larger mobile cations) to the improved cationic diffusion experimentally observed in honeycomb layered oxides. Furthermore, we elucidate the role of inter-layer distance and cationic size in cationic diffusion. Our theoretical studies reveal the dominance of inter-layer distance over cationic size, a crucial insight into the further performance enhancement of honeycomb layered oxides.
I. INTRODUCTION<!>II. METHODOLOGY A. Interatomic Potential<!>C. Estimation of key quantities<!>A. Framework Structure<!>B. Ionic Conductivity<!>C. Microscopic mechanisms of cationic transport<!>IV. CONCLUSION
<p>As the winds of change continue to push humanity to the crescendo of an energy revolution, high-energydensity storage systems have gained momentous traction, with lithium-ion batteries at the epicentre of both commercial and industrial applications. [1][2][3] This dominance of lithium-ion chemistries is principally propelled by their low redox potential (-3.04 V versus standard hydrogen electrode (SHE)) along with expedient ion kinetics that facilitate the development of layered oxide cathode materials that engender batteries with high voltages, highpower-densities and excellent cyclability. [2,3] However, the sustainability of lithium-ion batteries is heavily impeded by the prohibitive costs, safety issues and scarcity of lithium resources, galvanising explorations into alternative chemistries with terrestrial affluence such as sodium-ion and potassium-ion batteries. [4,5] Although the low redox potentials of Na and K (-2.71 V and -2.93 V versus SHE in aqueous electrolytes, respectively) earmark their potential as Li-ion alternatives, their layered oxide analogues have been found to yield low average voltages (∼ 3 V) due to their large ionic sizes. [6] In an effort to ameliorate the electrochemical performance and structural stability of layered oxides, partial substitution of the transition metal species with other transition metal species or highly valent species (such as chalcogens or pnictogens) has been employed. In this pursuit, a unique class of heterostructures known as honeycomb layered oxides has garnered prominence as high-energy-density cathode materials for their exquisite crystal structures that engender remarkable ion mobility and high voltages even with Na and K chemistries. [7][8][9][10][11][12][13][14][15][16][17][18][19] These materials typically adopt chemical compositions of 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 or coinage-metal atoms such as Cr, Mn, Fe, Co, Ni, Cu or some combination thereof; D represents pentavalent or hexavalent metal atoms such as Sb, Te, Bi; and A can be alkali atoms such as Li, Na, K, etc. [7,20] Given the differences in their ionic sizes and valency states, the 2:1 atomic ratio of M and D (such as in A 2 M 2 DO 6 and A 3 M 2 DO 6 ), renders 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.</p><p>Fast ion kinetics and high voltage capabilities are key prerequisites of any energy material geared towards future capacious battery applications. In this quest, the fine-tuning of the honeycomb layered structures which involves the judicious selection of resident species has been the focal point of numerous studies on this subject. [21][22][23][24][25][26] For high capacity, trivalent or divalent metals with high structural stability and high redox voltages (such as Ni 2+ /Ni 3+ /Ni 4+ in A 2 Ni 2 DO 6 , A 3 Ni 2 DO 6 or A 4 NiDO 6 ) are favoured as they allow partial alkali atom occupation without disintegrating. [8,10,[12][13][14][15][16][17][18][19]27] For high voltages, highly valent pnictogen or chalcogen atoms such as Te, Bi and Sb are employed to lower the covalency of the bond between the oxygen atoms and (di or tri) valent transition metals thereby increasing the energy required for (di or tri) valent cation oxidation. In turn, this induces a large increase in the overall voltage of the battery which can be predicted by the Pauling scale (Te > Sb > Bi) of electronegativity. Indeed, tellurates used alongside Ni 2+ in honeycomb layered oxides have been reported to produce the highest voltages (over 4V) to date. [8,10,12,13,[15][16][17][18] Thus, for a comparative study on honeycomb layered oxide electrochemical performance, Ni 2+ -based honeycomb tellurates (A 2 Ni 2 TeO 6 ) present an ideal pedagogical platform as a leverage for insights into the optimisation of related honeycomb layered oxides.</p><p>To surmise the mechanisms of honeycomb layered frameworks, a crystal structural illustration of A 2 Ni 2 TeO 6 (A = Li, Na and K) is furnished in Figure 1. As shown, the A atoms are interspersed between metal slabs consisting of Ni and Te in an octahedral coordination with oxygen atoms. The inter-slab (inter-layer) distances increase in the sequence of Li < Na < K in accord with the Shannon-Prewitt ionic size of the alkali atoms (Note that an ordered Li 2 Ni 2 TeO 6 polymorph [12] is used for ease of comparison with the Na and K counterparts). [15,21] As affirmed by previous structural studies on honeycomb layered frameworks, cations with smaller ionic radii (such as Li + in Li 2 Ni 2 TeO 6 ) are characterised by alkali planes with smaller inter-layer distances (consequently stronger inter-layer bonds) which require more energy to facilitate two-dimensional (2D) cation diffusion during (de)insertion processes. [7] On the other hand, resident cations with larger ionic radii (such as Na + and K + in Na 2 Ni 2 TeO 6 and K 2 Ni 2 TeO 6 ) foster weaker inter-layer bonds that facilitate facile 2D cationic diffusion during battery operations As such, Ni-based honeycomb tellurate compositions comprising Na and K have demonstrated faster cationic mobilities relative to their Li-based counterparts, with K displaying the lowest energy barrier for cationic mobility (viz., activation energy). [28,29] Despite the immense potential envisioned in the experimental and theoretical explorations of honeycomb layered materials, details surrounding the various mesoscopic mechanisms engendering their remarkable electrochemical performance still remain unexplored. Specifically, bigger size cations show higher activation energy barrier, resulting in lower cationic diffusion. This series of materials, A 2 Ni 2 TeO 6 (A = Li, Na and K), show completely contrasting behaviour. Thus, it attracts theoretical attention to understand the origin of such contrasting behaviour. In light of this, computational methodologies such as molecular dynamics (MD) simulations present a unique avenue to unveil the microscopic physicochemical properties such as ion transport and energetics at atomic scale, otherwise presently beyond experimental reach. [30] In fact, MD simulations have been shown to accurately depict the structural and transport properties of previously reported energy materials such as our previous works on the physicochemical properties of Nabased honeycomb layered tellurates (i.e., Na 2 M 2 TeO 6 , where M = Ni, Zn, Co and Mg). [31][32][33] Moreover, these simulations have been instrumental in garnering insights into the diffusive properties of other materials [34] that have later on been experimentally confirmed (e.g., in Na 2 Ni 2 TeO 6 ). [27,29] To accurately depict microscopic physicochemical properties of materials that are beyond experimental reach, it is therefore crucial to use reliable potential parameters that can effectively reproduce the structural and transport properties availed by experiment.</p><p>Thus, in this study we present an extensive MD study systematically performed using a refined set of interatomic potential parameters in A 2 Ni 2 TeO 6 (where A = Li, Na, or K). We focused on the fully discharged A 2 Ni 2 TeO 6 phases, since the crystallographic experimental data for the charged states of A 2 Ni 2 TeO 6 (i.e., A x Ni 2 TeO 6 (x < 2)) are lacking in literature. Nonetheless, distinct from most layered oxides, it is known that A 2 Ni 2 TeO 6 presents a finite number of vacancies, even in the fully discharged state (x = 2). [7,10,15,18,21,28] The MD simulations demonstrate the effectiveness of the Vashishta-Rahman form of the interatomic potential, [35] reproducing various structural and transport properties of this class of materials. Further, we predict that an increase in the inter-layer distance results in improved microscopic transport of cations through diffusion. As such, the larger ionic radii of K and Na relative to Li engenders an exponential enhancement of the cationic conduction of K 2 Ni 2 TeO 6 and Na 2 Ni 2 TeO 6 , in comparison with Li 2 Ni 2 TeO 6 . Through these investigations, we also connect the wider bottleneck radius formed by oxygens in cationic migration path for the larger mobile cations within the frameworks to the improved cationic diffusion experimentally observed in honeycomb layered oxides.</p><!><p>The previously reported Vashishta-Rahman form of the interatomic potential [35] is employed alongside parameters to produce the structural and dynamical properties of a series of honeycomb layered oxides, namely Na 2 M 2 TeO 6 (M = Ni, Zn, Co and Mg), [31] as follows; [35]</p><p>where q i represents the charge and σi ,σ j the ionic radius of the i-th, j-th ion. 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 ion pairs i and j, respectively. Compared to the more popular Born-Mayer (Buckingham) [36] and Lennard-Jones potentials, [37] 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) in particular, between the anion pairs. Other related parameters, not available in literature, were determined using empirical fitting to attain the experimentally reported bond lengths (Ni-O, Te-O, and O-O) of A 2 Ni 2 TeO 6 , where A = Li, Na, or K (Supplemental Material, Table SI). After sta-bilising the framework structure, the A-A and A-O parameters were refined to reproduce the experimentally reported cationic hopping activation energy, where A and O play the roles of the cations and framework anions respectively. [31,34,38,39] B. Computational Details</p><p>Employing the optimised set of parameters listed in Table S1, a series of MD simulations were carried out in the temperature range of 500 K to 700 K with constant temperature maintained at 25 K intervals, and zero atmospheric pressure. The Parrinello-Rahman isobaricisothermal (NPT) MD method, [40] which allows for changes in the simulation box sizes whilst keeping angles fixed was used. The temperatures and pressure in the system were controlled using thermostatting and barostatting techniques whereby some dynamic variables are coupled with particle velocities and simulation box dimensions. Simulations commenced from the structures derived from the respective X-ray diffraction patterns wherein the Li 2 Ni 2 TeO 6 simulation supercell was characterised by 3 × 5 × 2-unit cells comprising 1320 atoms in a orthorhombic symmetry and the A 2 Ni 2 TeO 6 (A = Na, and K)-systems consisted of 1100 atoms in 5 × 5 × 2 rhombohedral unit cells. For A = Na and K, several partially occupied crystallographic sites were identified. In the MD simulations, the respective Na and K ions were placed at the A2 sites of the starting structure to avoid strong cation-cation repulsion. For clarity to readers, the crystallographically distinct sites (A1, A2 and A3) in the layered frameworks of Na 2 Ni 2 TeO 6 and K 2 Ni 2 TeO 6 are shown Figure 1 . A time step of 1 fs was applied to the Velocity Verlet algorithm to solve the Newton's 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. In order to guarantee the thermodynamic convergence properties and ensure that the system size is adequate, a few longer run-time simulations of 100 ns as well as simulations on larger systems consisting of 6×6×3-unit cells were additionally performed, respectively. 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. [41] A cut-off distance of 10 Å was used for both the shortrange interactions and the short-range part of the Ewald summation. Canonical MD (NVT-MD) simulations were further performed using the average cell parameters obtained from NPT-MD, as NPT-MD modifies the kinetics to satisfy constant temperature and pressure conditions by sacrificing rigorous kinetic evolution. Nonetheless, we also performed microcanonical (NVE)-MD simulations for K 2 Ni 2 TeO 6 and compared it with our NVT-MD simulation results.</p><p>Additionally, we performed two sets of MD simulation corresponding to the following two scenarios: (i) a se-ries of NVT-MD simulations performed at 600 K with progressively increased inter-slab distances (only c-axis parameters) along with constant a, b -axes parameters, and ionic radius in order to understand the influence of inter-layer distance decoupled from ionic radius. (ii) A series of NVT-MD simulations conducted at 600 K with progressively increased cationic sizes, namely σ (1.20 Å to 1.30 Å at 0.05 Å intervals), whilst keeping the interlayer distance fixed to correspond with the ionic radius of K (i.e., 1.38 Å ), in order to understand the role of bottleneck radius formed by surrounded oxygen atoms. The aforementioned assumptions significantly divert from the already experimentally-verified correlation of inter-layer distance with ionic radius. Nonetheless, we should mention here that, for instance, subjecting this class of honeycomb layered materials to (positive) pressure corresponds to fixing the ionic radius whilst varying the inter-layer distance in accordance with scenario (i) above, whereas the converse (fixing the inter-layer distance and varying the ionic radius), in accordance with scenario (ii) above corresponds to expanding the material (i.e., inducing negative pressure). Consequently, conventional experimental results with unpressurised materials correspond to the zero-pressure scenario, where the inter-layer distance is strictly correlated with the ionic radius of the cations. As far as we can tell, experiments with finite pressure parameters have not yet been explored in experiments and simulations of this class of honeycomb layered oxides despite being standard in condensed matter experimental physics such as high-temperature superconducting materials where pressure modifies the transition temperature of the material and hence the conductivity, amongst other novelties. [42] In the case of the A 2 Ni 2 TeO 6 , pressure effects pertinent to the increase of inter-layer distance at fixed cationic size could be experimentally explored through chemical means, via the equation of state for the layered material. For instance, in a separate work, [43] two of the authors in this paper have conjectured that the curvature of the ab plane in these materials correlates with the number of cations/vacancies and/or the topology of the surface, which in turn requires the curvature to be proportional to the energy density (or more succinctly the number per unit volume) of the diffusing cations. Typically, the number per unit volume, ρ at thermodynamic equilibrium is related to pressure, P by an equation of state. [44] This avails the avenue to chemically change ρ or the curvature and hence P in the material. A more rigorous discussion has been appended in the Supplemental Material, section IB. In the simulation, we additionally point out that we fixed the inter-layer distance and separately the ionic radius to that of the K system as it possesses the widest inter-layer distance and the largest ionic radius amongst the honeycomb layered nickel tellurates.</p><!><p>Following Einstein's relation, the self-diffusion coefficient of Na + is calculated from mean square displacement (MSD) for the two-dimensional transport as,</p><p>where N t denotes the total number of atoms in the system, r j (t) is the position vector of the j th ion at time t and the angular bracket indicates average over several time origins. Note that the MSD trends against time remain unchanged for both NVT and NVE simulations, as is shown in Figure S1. The diffusion coefficient, D, depends on temperature (T ) following Arrhenius equation,</p><p>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 potential energy of individual cations is calculated as,</p><p>where V ij is the interaction potential given in eq. ( 1), such that the total potential energy of the system is,</p><p>Finally, the relative free energy (as seen by 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 by the formula (Supplemental Material, section IA),</p><p>where p represents the population density function (occupancy) of the cations, p max is the maximum value of the population density function, ∆V is the difference in potential energy of individual Also note that the population densities of cations remain unchanged for both NVE-and NVT-MD cases (Figure S2).</p><!><p>As shown in Figure 1, the A 2 Ni 2 TeO 6 compounds, where A = Li, Na, and K, have a layered structure of A cations sandwiched between transition metal  [15,21] whilst Li 2 Ni 2 TeO 6 is known to entail an orthorhombic and monoclinic crystal symmetries contingent on the synthetic protocols. [12] Moreover, all the octahedral layers in the Na 2 Ni 2 TeO 6 and K 2 Ni 2 TeO 6 are seen to be identical, whereas a slight shift between the top and bottom layer around the conduction plane is observed in the case of Li 2 Ni 2 TeO 6 (as shown in Figures 1(a), 1(b) and 1(c)). The layers are held together along the c-axis by Van der Waals forces and by interactions mediated via A atoms occupying their inter-layers. The refined set of interatomic potential parameters retains these structural features in conformity with experimental findings. The average cell parameters of A 2 Ni 2 TeO 6 calculated from NPT-MD simulations at 300 K are placed alongside the room-temperature experimental values in Table I, where the calculated caxis parameter (given the crystallographic arrangement of the layers in the c-axis direction) increases upon substitution of the A atom species in the order of Li, Na and K, respectively, in corroboration with experimental observations. [12,15,21] The progressive increase noted in the cell parameters of this family of materials is ascribed to increases in cation sizes in the order of Li, Na and K as a result of intensified ionic repulsion inside the layers as prescribed in the presence of larger sized cations. It should be noted that the deviations from experimental values in all cases are well within 3%.</p><p>For a comparison of the local structures, radial distribution functions, g(r), between the framework ion pairs, Ni-O, Te-O, and O-O, obtained from NPT-MD simulations at 300 K for the A 2 Ni 2 TeO 6 compounds (A = Li, Na, and K) are furnished in Figure 2 . For clarity, the radial distribution functions calculated from the respec-tive X-ray diffraction (XRD) structures [12,15,21,28] are also produced in the form of bar-plots. Noteworthy is that Li 2 Ni 2 TeO 6 shows structural disorders, resulting in the emergence of several XRD peaks compared to Na 2 Ni 2 TeO 6 and K 2 Ni 2 TeO 6 . Except for the thermal broadening, all the calculated peaks from MD simulations are reasonably consistent with their corresponding g(r), calculated for the XRD structure.For Te and Ni, exactly six oxygen coordination numbers are found within the 1st neighbouring distance, forming NiO 6 and TeO 6 octahedra. The average Te-O and Ni-O distances of the octahedra for A = Li, Na, and K, calculated from MD simulations (listed in Table II), also show good agreement with the reported XRD structure. Similarly, the radial distribution functions obtained from the other simulations performed under different conditions (see the METHODOLOGY section II) revealed that in all cases, the respective crystal structures were retained, affirming that the refined set of interatomic parameters successfully reproduce the structure.</p><!><p>The A cations are located in the ab-plane between the polyhedral layers and with loose coordinations to the framework layers which constitute a large number of cationic sublattices, specifically for A = Na, and K. These sublattices facilitate high ionic diffusion, in accordance with eq. (2a) and ascertained by the mean square displacements (MSD) of the cations plotted against time (Figure 3(a)). On the other hand, the closely packed octahedral layers parallel to the ab-plane restrict the cationic motion along the c-axis direction of the cell, as reflected in the inset of Figure 3(a) (less than 0.2 Å [2]). Therefore, it can be deduced that cationic mobility is restricted within the sublattices oriented parallel to the ab-plane, rendering the diffusion to be highly anisotropic (viz., two-dimensional). The MSD of the A ions is noted to increase drastically with increase in cationic size in the order of Li, Na and K. Figure 3(b) displays the diffusion coefficient, which was calculated from the slope of Figure 3(a) using eq. (2a), as a function of inter-layer distance. It is worth noting that the diffusion coefficient unequivocally follows the exponential dependency of the inter-layer distance, whereas the experimental activation energy diminishes linearly with increasing inter-layer dis-  tances as is displayed in Figure 3(b). [15,21] The logarithmic plots of D versus inverse temperature (1000/T ) are displayed in Figure 4(a). According to eq. ( 3), ln D = ln D 0 − E a /k B T corresponds to the equation of a straight line, where −E a /k B is the slope and ln D 0 the y-intercept. This implies that both the activation energy and the pre-exponential factor can be readily extracted from Figure 4(a). The activation energies were calculated using the diffusion constant (related to the diffusion coefficient in eq. ( 3)) instead of conductivity (Green-Kubo relation), an experimentally measurable quantity. It is worthy to note that conductivity, σ at thermal equilibrium and diffusion constant, D can be succinctly linked by σ = qµρ and µ = D/k B T , where µ is mobility, ρ is charge density and q is the unit charge. This implies that both approaches are equivalent in particular regimes (Supplemental Material, section IC). However, even though the Green-Kubo formula appears the more straightforward way to make the connection with measurable quantities such as conductivity, simulations with it are known to yield infinitely slow decay (or what is referred to as 'long-time tail') in conserved hydrodynamic variables. [45] This includes large fluctuations in the calculated electrical conductivity at long times. [46][47][48] The ) [21,28,43] and computation (MD simulations) against the inter-layer distances in A2Ni 2 TeO 6 (A = Li, Na and K). Note that the correlation plots of experimental against the computed ionic conductivity and activation energy values have been listed, demonstrating univocally that the Vashishta-Rahman potential adopted for our MD computation captures, with high fidelity, the appropriate physics of ion diffusion in this class of honeycomb layered oxides activation energy of Li 2 Ni 2 TeO 6 is found to be highest (0.56 eV) and lowest (0.26 eV) for K 2 Ni 2 TeO 6 amongst the three systems. A comparison of the activation energies calculated from the present MD simulations and reported experimental results is shown in Figure 4(b). Note that the MD results overestimate the activation energies for the three compounds compared to the reported ones. Usually, real material systems consist of several vacancies due to defects that lead to higher diffusion, which have not been considered here. Moreover, we have used the initial potential model of Na Alkali-ion diffusion in materials, in particular, entail the concomitant interplay between the inter-layer (interslab) distance and cationic size. In principle, a wider inter-slab distance reduces the interaction between the mobile cation and the ions in the transition metal slab, culminating in higher cationic diffusion (as identified in Figure 3(b), by the increase of cationic diffusion (diffusivity) coefficient (D) with increase in inter-layer distance. Consequently, due to the correlation of inter-layer distance with ionic radius in unpressurised materials, the cationic diffusion (diffusivity) is expected to increase with increasing cationic size. Therefore, it is difficult to ascertain the role of cationic size since the enhancement of diffusivity in K systems can also be attributed to the increased inter-layer distance. Fortunately, it is possible in theoretical simulations to decouple the inter-layer distance from the cation size and treat them independently of each other even though they are typically correlated in unpressurised systems. Therefore, in the MD simulations, the c-axis parameter, which corresponds to the inter-layer distance, is increased whilst the cation size is kept constant. The c-axis parameter was progressively increased from 11.3 to 13.3 Å at intervals of 0.2 Å , to cover the inter-layer distance range from Li to K. The ionic radius was kept constant at the value corresponding to that of K (i.e., 1.38 Å ), as K-ions have the largest ionic radius reported so far for this family of honeycomb layered oxides. As can be seen in the plot of the inter-slab distribution (furnished as Figure S3 in the Supplementary Information), increasing the c-axis parameter basically increases the inter-slab distance, d (i.e., d = c/2) without forming other structural artifacts in the MD simulation, such as undulation. However, it is important to note that different size cations engender different interlayer distances to maintain structural stability. Thus, through MD simulations, it is possible to study microscopic phenomena beyond the reach of practical experiments.</p><p>Accordingly, the inter-layer distances below the ionic radius of lithium were not considered because lithium cations are the smallest cations that embody the present honeycomb layered framework. Moreover, the cationic diffusion was observed to diminish with decrease in the inter-layer distance (c-axis parameter), as shown in Figure S4 ; indicative that inter-layer distances smaller than 11.3 Å , would result in unfavourable cationic diffusion. Notably, the diffusivity trend significantly deviates from exponential dependency of inter-layer distances as shown in Figure 5, whereby the logarithmic of diffusion coefficient plot shows a converging nature at higher inter-layer distances. Although this might be a consequence of the current simulation conditions, where the a-and b-axes parameters and the ionic sizes are kept constant, it might also be an indicator that larger inter-layer distances do not necessarily result in higher cationic diffusion presumably as a result of structural instability or diminished inter-layer interactions.</p><!><p>Although the aforementioned results are already illuminating, it is imperative to take into consideration the atomistic mechanisms governing the cations confined within the sublattices along the ab-plane in order to accurately discern the nature of the enhanced diffusion in Na 2 Ni 2 TeO 6 and K 2 Ni 2 TeO 6 compared to Li 2 Ni 2 TeO 6 . On that account, three different cationic sites were identified from experimental studies inside the conduction plane of A 2 Ni 2 TeO 6 (A = Na and K) in each inter-layer of a unit cell, labelled as A1, A2 and A3, with a multiplicity of 3, 2, and 1 respectively. [15,21,27,28] (For clarity to readers, A1, A2 and A3 are crystallographically distinct sites.) These sites are the centre of trigonal prisms formed by oxygen atoms, as shown in Figure 1. The octahedral layers include tetrahedral voids (or 'holes') formed between two neighbouring NiO 6 and TeO 6 octahedra with the A1 site located mid-way between the two tetrahedral voids from the top and bottom layers. The A2 sites are sandwiched between the triangular faces of NiO 6 octahedra from the top and bottom layers, whereas the A3 site is located between two TeO 6 octahedra. However, the experimentally reported Li 2 Ni 2 TeO 6 shows a different structure whereby the NiO 6 octahedra are slightly deformed, and the layered stacking sequence is different from the Na 2 Ni 2 TeO 6 and K 2 Ni 2 TeO 6 which have the TeO 6 octahedra at the top and bottom of the A3-sites. The smaller inter-layer separation (5.08 Å ) in Li 2 Ni 2 TeO 6 (compared to Na 2 Ni 2 TeO 6 (5.57 Å ) and K 2 Ni 2 TeO 6 (6.23 Å )), engenders a high Te-Te Coulombic repulsion along the c-axis direction which prevents the TeO 6 octahedra from achieving a stable vertical orientation. Thus, unlike the Na and K systems, the Li cations orient themselves at the top and bottom of the octahedral voids in order to avoid the resultant strong repulsion at the A3 site.</p><p>In order to understand cationic site preferences and hopping energy barriers between sites, the population density (number per unit area) pattern and average potential energy of the mobile cations were calculated in accordance with eq. ( 5) (detailed in the METHODOL-OGY section II). Both the population density and potential energy in the conduction plane are projected onto a 2D grid on the ab-plane (with all the coordinates of the cations folded into a single cell). Figure 6 displays the population density, whilst Figure 7 shows the corresponding potential energy surface of cations in the order of Li, Na and K (replicated over 2×2-unit cells for easy visualisation of their continuity). The Li population densities are mostly localised due to the deep potential wells with minima of -3.2 eV (Figures 6(a) and 7(a)), whereas the Na-substituted system exhibits well-connected cationic diffusion channels between the A1 and A2 sites (Figures 6(b) and 7(b)). Due to the aforementioned strong repulsion, there are no Na cations occupying the A3 sites, which are located directly between the top and bottom Te-Te octahedra along the c-axis direction. In the case of Na 2 Ni 2 TeO 6 , the potential wells are relatively shallow (-2.6 eV) compared to the Li system. Conversely, K 2 Ni 2 TeO 6 exhibits slightly different profiles (Figures 6(c) and 7(c)), characterised by the shallowest potential wells (-2.2 eV) amongst the three systems. Thus,  6 are related to the potential energy profiles via eq. ( 6). Thus, the difference in the two profiles can be attributed to a finite relative entropy, ∆S in the system. even though the population density profile of cations in K 2 Ni 2 TeO 6 is slightly different than of Na 2 Ni 2 TeO 6 , the path connectivity of K + is essentially the same as that simulated for Na + , indicating the same cation hopping route. The shallower potential profile in the K system attests to its considerably lower potential energy barrier, allowing for more facile shifts of K cations to and from the A1 and A2 sites compared to the Na system. The calculated potential energy difference from A1 to A2 is 0.25 eV; albeit a slightly different value than reported based on bond valency energy landscape (BVEL) summation. [18] Note that BVEL heavily relies on experimental data for tuning, thus compromising its accuracy and predictive power. Thus, the results in Figure 7 evince that the minima values of the potential energy increase (i.e., the potential becomes shallower) in the order of Li, Na and K, which also corresponds to the increasing order of the inter-layer distances, indicative that the potential energy barrier inhibiting the hopping of the cations be-tween sites reduces with inter-layer distance, leading to a drastically higher cationic diffusion. This was additionally affirmed by calculations of the profiles in relation to the progressively increasing inter-layer distances, as shown in Figure 5.</p><p>For a more detailed quantitative estimation of the role of the cationic occupancy and its associated (relative) free energy landscape, we considered the 2D population density distribution for progressively increasing c-axes parameters (Figure S5), and calculated the (relative) free energy (Figure S6(b)) from the population distribution (Figure S6(a)) along a straight-line path conjured up by an imaginary cylindrical axis with a radius of 0.8 Å connecting the neighbouring A2 and A1 sites in accordance with eq. ( 6) (see the METHODOLOGY section II). Note that, the Na and K systems exhibit slightly differing A1 and A2 site occupancy, as is revealed by the population density profile shown in Figures 6(b) and 6(c) respectively. For instance, in the Na system, the population densities at the A1 and A2-sites are almost identical whereas in the K system, the population at the A2-sites is significantly higher than the A1-sites. Remarkably, (relative) free energy distribution calculations from their respective population density distributions revealed that although K system exhibits a lower activation energy than the Na system, the (relative) free energy barrier for K ion diffusion (0.07 eV) is higher than its Na counterpart (0.04 eV)(Figure S6(b)). This juxtaposed behaviour can be ascribed to the site topology and the connectivity of the two systems. Moreover, as illustrated in Figure 6(b, c) and Figure 7(b, c), the A1 sites are connected by two A2 sites, whereas A2 sites are connected by three A1 sites implying that a cation occupying an A1 site has two possible paths (corresponding to a configuration entropy of ∼ k B ln 2), whereas it finds three possible paths if it occupies an A2-site (configuration entropy of ∼ k B ln 3). Therefore, because the Na system is characterised by higher population densities at the A1-sites, cationic diffusion from A2 to A1 will be restricted by the lower availability of paths (lower configuration entropy). Conversely, a lower population density implies a higher number of cationic vacancies. Hence the vacancies at A2 are larger than at A1, and by extension, Na cations can hop from A1 to A2, corresponding to a 'hole' hopping from A2 to A1 sites, where the 'holes' also have three available A2 to A1 pathways (configuration entropy of ∼ k B ln 3), and hence contributes an opposite diffusion current. Since the pathway with the larger configuration entropy contributes to the bulk of the diffusion current in the system, it implies that in the Na system, the 'holes' will contribute the bulk of the diffusion current, thus offsetting any low configuration entropy restrictions to the diffusion current imposed on the cations in the A1 sites. Likewise, the lower population density of K cations at A1 sites compared to A2 sites and the lower configuration entropy restricts the 'holes' from hopping from the A1 sites to the A2 sites whilst the higher population density and configuration entropy of the K cations at the A2 compared to the A1 site avails more avenues for the K cations to hop from A2 sites to A1 sites, implying that the bulk of the current in K systems is as a result of the high configuration entropy of the cations at A2 site. Thus, this cation-'hole' picture also implies that the cation occupancy/vacancies and population density distributions alone (Figure 6(b) and Figure 6(c)) cannot explicitly explain why the K system exhibits, in experiments, better diffusion properties than the Na system. In fact, this places more emphasis on the role of activation energies (Figure 4(b)) and hence their diffusion coefficients, as well as the difference in the bottleneck radius (Figure 8) of the two systems, in order to ascertain this trend in experiment.</p><p>For a deeper understanding of the role of the interlayer distance in cationic diffusion at an atomistic level, the bottleneck radius formed by the oxygen atoms, was calculated along the cationic diffusion pathway. As defined by Evstigneeva and co-workers, [21] the bottleneck FIG. 8. Schematic representation of the cation-hopping path from A2 to A1 sites across the A2Ni 2 TeO 6 topology and the associated (relative) free energy landscape. The bottleneck radius and the cationic migration path have been highlighted. A curve illustrating the free-energy landscape is also appended, for clarity. The quadrangular bottleneck is located near the transition state (one of the bottleneck radii is shaded for emphasis). Different colours have been used for atoms residing in A1 and A2 sites, for clarity.</p><p>radius is the smallest A-O distance at the saddle point in the cationic diffusion path from the A1 to A2 sites, as shown in Figure 8. In the present study, the other bottlenecks formed between the path A1 to A3 are not taken into account, as they do not play any role on cationic diffusion owing to the cationic absence at the A3-sites. The results indicate that the bottleneck distance for the Na system (2.35 Å ) is smaller than the K system (2.57 Å ) as is evident in Figure 9 . Thus, the Na system embodies a stronger force of attraction at the saddle points of the (relative) free energy path owing to the closer O atoms, effectively inhibiting cationic diffusion along the channels. Conversely, the K system manifests the weakest force of attraction at the saddle points due to the distant O atoms, yielding the highest cationic diffusion rate as observed in the MD simulation results. Figure 9 indeed depicts a similar trend. That is, larger inter-layer distances confer wider bottleneck radii; thus, resulting in higher cationic diffusion, as further displayed in Figure 5. By varying ionic radius whilst keeping the inter-layer distance fixed to 6.23 Å , corresponding to the value reported in measurements of the unpressurised K 2 Ni 2 TeO 6 , the effective bottleneck radius increases the smaller the cationic size. In other words, for a fixed bottleneck radius, a smaller-sized cation can easily pass through (compared to a larger-sized cation), thus resulting in higher cationic diffusion (as is evident in Figure S7). This behaviour is known as 'levitation effect', as studied elsewhere. [49] However, the diffusion FIG. 9. Bottleneck radius formed by oxygen atoms located between the cationic path from A1 to A2 for the different c-axis parameters. The bottleneck radii for Na and K are marked using blue and green square symbols, respectively. behaviour does not change as significantly as in the case of the change in inter-layer distance. Therefore, this clearly indicates the domination of inter-layer distance over the cationic size in establishing the leading order trends in the diffusion behaviour. The accurate depiction of the progressive increase in the bottleneck radius with increasing inter-layer distances reaffirms the use of MD simulations as an auspicious tool to decipher the atomistic mechanisms of honeycomb layered oxides.</p><p>In light of the observations presented in this study, it can be deduced that the de novo design of chemical compositions of A 2 Ni 2 TeO 6 entailing lithophile atoms with larger ionic sizes than K such as A = Rb or Cs could be a propitious strategy for enhancing the cationic diffusion for honeycomb layered oxide family discussed herein. However, MD simulations suggest that such systems are relatively unstable, because of the extremely weak Van der Waal interactions between adjacent layers. Stabilisation of such new compositions may necessitate the utilisation of high-pressure synthesis. Moreover, computation has also aided us to understand the high K-ion diffusion experimentally reported in K 2 Ni 2 TeO 6 using muon spin rotation and relaxation spectroscopy. [28] The high K-ion diffusion is explained from the potential energy barrier (weaker cation and inter-layer interaction), favourable path connectivity and wider bottleneck radius.</p><!><p>This study shows the efficacy of using the Vashishta-Rahman potential in reproducing structural and transport properties of honeycomb layered oxides based on nickel tellurates and this potential model can be leveraged to garner insights on related honeycomb layered oxide materials. Moreover, the ability to depict, with high fidelity, the atomic scale cationic diffusion aspects of honeycomb layered oxides, avails a new avenue to test the ion dynamics of various honeycomb layered oxides when subjected to various conditions (e.g. pressure) presently unreported and currently inaccessible via experiment.</p><p>The interplay between inter-layer distance and cationic size in the ionic diffusion of honeycomb layered oxides is unveiled -a wider inter-layer distance favours faster diffusion, whereas a larger cationic size does not necessarily favour faster diffusion. Based on the computational results, cationic diffusion is influenced more by the inter-layer distance than the size of the cation. Moreover, we note that, there is an exponential increase in the ionic diffusion with increasing inter-layer distance, which puts forward the possibility to test various theoretical models [43] relating to ionic diffusion in honeycomb layered oxides as well as stimulate new experimental avenues that decouple the effects of inter-layer distance from the ionic radius such as applying pressure or recreating an expansion scenario to test the predicted trends herein for negative pressure.</p><p>Therefore, the results reported herein not only affirm the remarkable performance of honeycomb layered oxides but also validate the use of computation techniques such as molecular dynamic simulations as a propitious path in the quest to design high-performance materials for future capacious batteries.</p><p>V. ACKNOWLEDGEMENTS</p>