Exploring superionic conduction in lithium oxyhalide solid electrolytes considering composition and structural factors

Crystal structureThe reported series of Li2+xZrCl6-xOx (LZCO) materials were mechanochemically synthesized with multiple phases depending on oxygen content. The crystal structures of Li2+xZrCl6-xOx are composed of different sites, intrinsic vacancies, and partial occupancies, which makes these structures complex for theoretical study. LZCO with both \({\rm{P}}\bar{3}{\rm{m}}1\) and C2/m phases were considered for the simulations. The \({\rm{P}}\bar{3}{\rm{m}}1\) phase contains two different Li atoms, Li1 and Li2 occupying 6 g and 6 h sites, four different Zr atoms, Zr1, Zr2, Zr3 and Zr4 at 1a, 2d, 1b and 2d sites, three Cl and O atoms at three different 6i sites with partial occupancies, respectively. C2/m phase consists of three different Li atoms Li1, Li2 and Li3 occupying 4 h, 2d and 4 g sites, one Zr atom at 2a sites while two Cl and O sharing sites at 8j and 4i sites respectively. The occupancy of each atom at each site is different and directly related with O content in every composition.The phases \({\rm{P}}\bar{3}{\rm{m}}1\) and C2/m of Li2+xZrCl6-xOx (x = 0, 0.25, 0.5, 0.75, 1, 1.25,1.50 and 2) were simulated to extract their reflection of powder diffraction plotted in Fig. 1. With x = 0, Li2ZrCl6 for the phase \({\rm{P}}\bar{3}{\rm{m}}1\) has the highest peak near 30° of 2θ and the second highest peak near 15° of 2θ, however the first highest peak starts to increase near 15° of 2θ while second highest peak is near 30° of 2θ° for x = 0.25. The simulated pattern for 0.5 ≤ x ≤ 2 similar pattern was observed with the first highest peak near 15° of 2θ° and the second highest peak near 30° of 2θ°, whereas the 3rd and 4th peaks observed at 42° and 50° of 2θ° are same for all LZCO compositions (Fig. 1a). The pattern exhibits pure phase and isostructural with Li2ZrCl6 about the composition of O with x ≤ 0.25 and is in agreement with the experimental observed XRD pattern (Supplementary Fig. 1)31. The simulated pattern for the C2/m phase for 0 ≤ x ≤ 0.75 indicates no change in the peaks (Fig. 1b). The crystallographic data for x = 0.50 and x = 0.75 (Li2.5ZrCl5.5O0.5 and Li2.75ZrCl5.25O0.75) for C2/m phase provided in the experiment is inconsistent with the chemical formula, which may be the reason for no change in peaks. At x = 1, the first highest peak for C2/m phase is observed near 15° of 2θ° and the second highest peak near 30° of 2θ°, suggesting a phase transition. At x = 1.25 there is change and increase in the first highest peak and decrease in the second highest peak, which may be attributed to the coexistence of \({\rm{P}}\bar{3}{\rm{m}}1\) and C2/m phases. There is again change in the first and the second peaks for x = 1.5, and much decrease in the peaks between 10° and 20° of 2θ° and again increase in peak near 20° of 2θ° in C2/m phase (Fig. 1b). These simulated diffraction patterns indicate the discrepancy between the crystal structures in C2/m phase with the increasing oxygen doping content. The reflection of crystallinity and amorphous behavior for simulated XRD pattern of Li2+xZrCl6-xOx (x = 0, 1 and 2) are further analyzed. The selected structures indicated maximum crystallinity with negligible amorphous phase (Supplementary Fig. 2).Fig. 1: Simulated XRD patterns and electrostatically most stable structures.Simulated XRD patterns of Li2+xZrCl6-xOx for phases with a space group \({\rm{P}}\bar{3}{\rm{m}}1\) and b space group C2/m. Crystal structure of Li2.5ZrCl5.5O0.5 with \({\rm{P}}\bar{3}{\rm{m}}1\) phase c with actual occupation of all atoms at each site, d electrostatically most stable structure with experimentally reported occupation sites, e electrostatically most stable structure with deviation in occupation sites, and f atomic and volumetric relaxed structure after spin polarization calculation. C2/m phase g with actual occupation of all atoms at each site, h electrostatically most stable structure with experimentally reported occupation sites, i electrostatically most stable structure with deviation in occupation sites, and j atomic and volumetric relaxed structure after spin polarization calculation. Cl/O indicates O doping at Cl sites, Li/Vac indicates partial occupancies of Li ion at different sites with Li intrinsic vacancies.The structures of Li2+xZrCl6-xOx are complex with intrinsic vacancies or partial occupation at all sites. Energetically most stable structures are needed to compute the key properties of electrolyte like phase stability, thermodynamic stability, and ionic transport. To have detailed structural analyses, Li2.5ZrCl5.5O0.5 with \({\rm{P}}\bar{3}{\rm{m}}1\) and C2/m phases were considered for comprehensive study using experimental crystallographic information31. Supercell program35 is useful tool which helps to generate all possible structures based on their site occupancies using permutations and combinations. To maintain the crystal charge neutrality and minimum number of atoms meeting the reported chemical formula, supercell program was used to compute electrostatic Coulomb energy. Electrostatic Coulomb energy of each structure is computed using valency and site occupation of each atom (Li1+, Zr4+, Cl1-, and O2-). The experimentally reported Li2.5ZrCl5.5O0.5 \({\rm{P}}\bar{3}{\rm{m}}1\) phase have Li1, Li2, Cl1, Cl2, Cl3, O1, O2, O3, Zr1, Zr2, Zr3 and Zr4 with occupancy of 0.53, 0.62. 0.95, 0.923, 0.977, 0.05, 0.077, 0.023, 0.649, 0.565, 0.709 and 0.256, respectively, while C2/m phase having Li1, Li2, Li3, Cl1, Cl2, O1, O2, and Zr1 with occupancy of 0.58, 0.2. 0.77, 0.825, 0.90, 0.175, 0.10 and 1.00, respectively (Fig. 1c, g). A series of crystal structures were generated, and the electrostatically most stable (lowest Coulomb energy) structure were extracted (Fig. 1d, h) for both phases. To make crystal structure more reliable, a little deviation in sites for both phases was made (\({\rm{P}}\bar{3}{\rm{m}}1\) phase having Li1, Li2, Cl1, Cl2, Cl3, O1, O2, O3, Zr1, Zr2, Zr3 and Zr4 with occupancy of 0.50, 0.75. 0.9167, 0.8334, 1.00, 0.0834, 0.1667, 0.00, 0.50, 0.50, 0.50 and 0.50, respectively, while C2/m phase having Li1, Li2, Li3, Cl1, Cl2, O1, O2, and Zr1 with occupancy of 0.25, 0.50. 0.75, 0.825, 0.90, 0.125, 0.0 and 1.00, respectively) and electrostatically most stable structure were obtained (Fig. 1e, i). Ten electrostatically lower energy structures from each phase were considered to get energetically most stable for which spin polarization calculations were carried out to get atomically and volumetrically relaxed structures (Fig. 1f, j). Coulomb energy for each series and phase is plotted in Supplementary Fig. 3. The ground state energy structure is used to compute optimized lattice parameters for each structure. Supplementary Fig. 4a indicates that the structure number 53 generated by supercell program has lowest electrostatic Coulomb energy (−311.564 eV), lowest total energy (−127.347 eV), and less hull energy (33 meV) computed by DFT calculation (Supplementary Table 3). These initial calculations are performed on the unit cell structure (28 atoms) with lattice parameter a = b = 10.93700 Å, c = 6.02200 Å for P3m1 phase. The structure number 53, which has minimum ground state energy is used to compute optimal lattice parameters, Supplementary Fig. 4b indicates that the optimal volume is little lower than as compared to experimental reported lattice parameters. A similar process is used for C2/m phase for all model structures (including LiCl deficient structure) used in this study. A supercell of 1 × 1 × 2 with 56 atoms and 2 × 1 × 2 supercell containing 76 atoms of \({\rm{P}}\bar{3}{\rm{m}}1\) and C2/m phases respectively, whereas forty-eight (48 atoms) atomic model of 1 × 1 × 2 of \({\rm{P}}\bar{3}{\rm{m}}1\) for amorphous structure was used for computing diffusion coefficients and ionic conductivity.Phase stability and electrochemical windowThe spin polarized calculations were performed to compute the stability, hull energy and electrochemical stability of the pure Li2ZrCl6, Li2.5ZrCl5.5O0.5 and Li1.75ZrCl4.75O0.5. It was found that both phases (\({\rm{P}}\bar{3}{\rm{m}}1\) and C2/m) of Li2ZrCl6 are thermodynamically stable with zero hull energy (equation 1), which is reflected by the ternary phase diagram at 0 K (Fig. 2a). However, oxygen containing structures (Li2.5ZrCl5.5O0.5 and Li1.75ZrCl4.75O0.5) are thermodynamically unstable with positive hull energy at 0 K and can be decomposed into other equilibrium phases according to following reaction:$${{\rm{Li}}}_{2}{\rm{Zr}}{{\rm{Cl}}}_{6}\to {{\rm{Li}}}_{2}{\rm{Zr}}{{\rm{Cl}}}_{6}$$
(1)
$$4\,{{\rm{Li}}}_{2.5}{\rm{Zr}}{{\rm{Cl}}}_{5.5}{{\rm{O}}}_{0.5}\to 3{\rm{Zr}}{{\rm{Cl}}}_{4}+10\,{\rm{LiCl}}+{{\rm{ZrO}}}_{2}$$
(2)
$$4\,{{\rm{Li}}}_{1.75}{{\rm{ZrCl}}}_{4.75}{{\rm{O}}}_{0.5}\to 3\,{{\rm{ZrCl}}}_{4}+7\,{\rm{LiCl}}+{{\rm{ZrO}}}_{2}$$
(3)
Fig. 2: Phase diagrams, electrochemical stability window voltages and total density of states.Ternary and quarterly phase diagram constructed at 0 K for a Li2ZrCl6 and b Li2.5ZrCl5.5O0.5. Electrochemical stability window voltage for c Li2ZrCl6, and d Li2.5ZrCl5.5O0.5. Total Density of states computed for Li2.5ZrO5.5Cl0.5 using e non-spin-polarized calculation with GGA functional, f spin-polarized calculation with GGA functional, and g non-spin-polarized calculation with HSE06 functional.The structure of Li2.5ZrCl5.5O0.5 with \({\rm{P}}\bar{3}{\rm{m}}1\) is energetically more stable with hull energy of 33 meV/atom and 41 meV/atom for C2/m phase, however both phases with O are thermodynamically unstable. Figure 2b shows the quarterly phase diagram of Li2.5ZrCl5.5O0.5 with green square in the inset green triangle indicating that it can be decomposed into the stable phases of ZrCl4, LiCl and ZrO2. Moreover, Li-depletion is an effective strategy to introduce structural variation and produce lattice defects for enhanced ionic transport in inorganic solid electrolytes. Li1.75ZrCl4.75O0.5 is a defect structure with the intrinsic depletion of Li0.75Cl0.75 in the Li2.5ZrCl5.5O0.5. However, the intrinsic deficiency of Li0.75Cl0.75 makes the structure even more unstable with the hull energy of 74 and 90 meV/atom in the \({\rm{P}}\bar{3}{\rm{m}}1\) and C2/m phase, respectively. In general, such type of structures with hull energy more than 50 meV/atom are not easy to be synthesized but are accessible under certain circumstances36,37,38. It is mentioned above that 300 crystal structures were generated for Li2.5ZrCl5.5O0.5 in this study by using the experimentally reported occupation of each atom, while 600 structures were extracted by deviation in occupation which exhibit low hull energy. However, no crystallographic information has been provided in the experimental study for deficient Li1.75ZrCl4.75O0.5. The modified structure information for Li2.5ZrCl5.5O0.5 and Li1.75ZrCl4.75O0.5 is provided in Supplementary Tables 1 and 2 of Supplementary note 1.The electrochemical stability window against Li electrode for Li2ZrCl6, Li2.5ZrCl5.5O0.5 and Li1.75ZrCl4.75O0.5 was investigated using Li grand potential phase diagram. Li2ZrCl6 exhibits a wide electrochemical stability window voltage (1.66−4.29 V versus Li/Li+) as compared to Li2.5ZrCl5.5O0.5 (1.75−3.96 V) indicated by the green line in Fig. 2c, d, respectively. However, the results suggest that both Li2ZrCl6 and Li2.5ZrCl5.5O0.5 are not stable with Li metal anode. The reduction and oxidation potentials and decomposition reactions into stable products are given in the Supplementary note 2. The electrochemical stability window voltage of Li1.75ZrCl4.75O0.5 is same with that of Li2.5ZrCl5.5O0.5 (Supplementary Fig. 5) for \({\rm{P}}\bar{3}{\rm{m}}1\) and C2/m phases.The electronic property analysis was further carried out to compute bandgap for Li2.5ZrCl5.5O0.5 and Li1.75ZrCl4.75O0.5 using spin, non-spin polarized calculations with GGA functional and HSE06 functional. Supplementary Fig. 6b shows that Li2.5ZrCl5.5O0.5 is a large bandgap (≈5 eV) insulator that could be used for practical and dendrite free electrolytes for batteries39,40, however the peaks of EDOS are not smooth. To make the EDOS smooth, Gaussian smearing method is used to extract the EDOS for GGA and HSE06 method with average number of configuration of 1000 and smearing value of 0.1 eV. The bandgap values computed by GGA functional (Fig. 2e, f) are much lower than bandgap of HSE06 functional (Fig. 2g), which results in the underestimation of GGA functional. Total and partial density of states (Supplementary Figs. 6 and 7) were also examined for both phases, with all the compositions exhibiting wide bandgaps (~5 eV).Li-ion transport propertiesTo examine ionic transport, AIMD simulations were performed using NVT ensemble. Mean square displacement (MSD) and diffusion coefficients were obtained by the linear fitting of MSD from Vaspkit41 for Li2.5ZrCl5.5O0.5 and Li1.75ZrCl4.75O0.5. The Li2.5ZrCl5.5O0.5 structure exhibits high Li MSD along three directions for both phases (Fig. 3a, c), which indicates fast three-dimensional diffusion of Li. The activation energy and Li ionic conductivity were examined by the Arrhenius plot of diffusivity (Fig. 3b, d). \({\rm{P}}\bar{3}{\rm{m}}1\) phase of Li2.5ZrCl5.5O0.5 exhibits a high Li-ionic conductivity of 1.63 mS cm-1 with an activation energy of 0.29 eV, while C2/m phase exhibits a lower conductivity of 0.24 mS cm−1 with an higher activation energy of 0.36 eV, which is in agreement with experimentally reported ionic conductivity31. No diffusion of Zr, Cl and O was observed during AIMD simulations (Supplementary Fig. 8), suggesting that the Li2.5ZrCl5.5O0.5 structure maintains its crystallinity at higher temperatures as well, whereas deficient Li1.75ZrCl4.75O0.5 could not maintain its crystallinity at temperature as low as 600 K, and melted (Supplementary Fig. 9). These results show that Li1.75ZrCl4.75O0.5 is disordered structure with excess depletion of 0.75LiCl and makes it energetically unfavorable, which is investigated below in the amorphous part.Fig. 3: MSD and Arrhenius plots.Total and directional Li MSD extracted for Li2.5ZrCl5.5O0.5 at 900 K for a \({\rm{P}}\bar{3}{\rm{m}}1\) phase and b C2/m phase. Arrhenius plots of diffusivity at different temperatures (600−1000 K) to compute extrapolated Li ionic conductivity at room temperature for Li2.5ZrCl5.5O0.5 c \({\rm{P}}\bar{3}{\rm{m}}1\) phase, and d C2/m phase.Pair correlation function (PCF) g(r) of Li2.5ZrCl5.5O0.5 was extracted for more comprehensive and structural analysis from AIMD simulations and Li ion transport mechanism. The PCF between Li-Cl (Fig. 4a) and Li-O (Fig. 4b) is presented, which indicates that Li-O pair is the main contributor in highest peaks of PCF, suggesting, O is strongly arranged at their sites during Li diffusion. The first highest peak of Li-Cl is at the radius of 2.4 Å, while Li-O first peak is at 1.9 Å at each simulation temperature for 100 ps. The second highest peak for Li-Cl and Li-O is at the radius of 4.6 Å and 4.4 Å, respectively. The inset Fig. 4a, b correspond to the PCF results, where solid lines represent bond lengths associated with the first peaks, while dotted lines denote the second highest peaks and beyond. The height of the peaks for both pairs decrease with the increase in temperature from 600 K to 1000 K, which suggests enhanced Li diffusion with rise in simulation temperature. The PCF files between Li-Li, Li-Zr, Li-Cl, and Li-O for Li2.5ZrCl5.5O0.5-C2/m phase was also calculated and presented in Supplementary Fig. 10 at each temperature, which indicates similar behavior as observed in \({\rm{P}}\bar{3}{\rm{m}}1\) phase. The Li-O pair has the maximum peak and decreased with the increasing temperature (Supplementary Fig. 11), suggesting high Li diffusion with increasing temperature.Fig. 4: Pair correlation function and vibrational density of states.Pair correlation function extracted from AIMD simulations at each temperature (600−1000 K) for 100 ps of Li2.5ZrCl5.5O0.5 for a Li-Cl pair and b Li-O pair. Vibrational density of states (VDOS) obtained from AIMD simulations at each temperature (600−1000 K) for 100 ps of Li2.5ZrCl5.5O0.5 of c Li and d O atom.VDOS was obtained from AIMD simulations versus frequency (eV) at temperatures of 600−1000 K for Li (Fig. 4c) and O (Fig. 4d) in Li2.5ZrCl5.5O0.542,43. The VDOS of Li at each temperature shows the highest peak near frequency of 0.03 eV and becomes straight after 0.1 eV. There is a decrease in peak with increase in simulation temperature (Supplementary Fig. 12), which is again an indication of fast Li diffusion with rise in temperature. As a case study for immobile species, VDOS is computed for oxygen. It can be seen in Fig. 4d that the peak values are not normalized at a single frequency, indicating O is strongly bonded and resist the movement during AIMD simulations at each temperature. The result suggests that O remained localized at their mean position and no diffusion is observed for the immobile Zr, Cl and O, which confirms the MSD result during AIMD simulations.Amorphous structureIn order to determine the formation of amorphous phase and ionic conduction mechanism, the Li1.75ZrCl4.75O0.5 structure is simulated using \({\rm{P}}\bar{3}{\rm{m}}1\) phase as an example. The structure was initially heated at 100 K for short simulation time, and heated at 300 K for 50 ps, at 600 K for 100 ps and finally heated at 1000 K for 100 ps to completely melt the structure. The structure was optimized after 1000 K with spin polarized DFT calculation. The optimized structure was finally quenched to 100 K with a constant rate of temperature gap of 200 K (800, 600, 400, 200 K) for 10 ps at each temperature using a time step of 2 fs with NVT ensemble using Nose‘-Hoover thermostat44. We have adopted the same method for the formation of amorphous just like other halide amorphous structures reported, including MAlCl2.5O0.75 (M = Li, Na)45, NaTaCl627, LiTaCl646, antiperovskite Na3-xOHxCl47. The complete melting, cooling and fully equilibrated process is plotted in Fig. 5a. AIMD simulations were carried out to heat the structure at 1000 K for 80 ps to completely melt the structure. It is a fact the structure shows diffusion for Cl, Zr and O at 600 K as well, however 1000 K is chosen so that the bonding of atoms, position of all atoms completely changes and melts. Finally, the structure was heated at different temperatures like 340, 360 and 380 K to get extrapolated ionic conductivity at 300 K. The structure was also directly simulated at 300 K to get direct Li diffusion and ionic conductivity. MSD and diffusion coefficients were extracted at 340, 360 and 380 K for 80 ps simulations time for Li1.75ZrCl4.75O0.5 (Fig. 5b). Extrapolated room temperature ionic conductivity (43.3 ± 3.3 mS cm−1) was obtained by Arrhenius plot (Fig. 5c) with an activation energy of 0.25 ± 0.1 eV. During heating at 1000 K the structure was completely melted with a high rate of diffusion of each element (e.g., Cl MSD plotted in Supplementary Fig. 13a). Although amorphous phase is obtained after cooling at 1−300 K, however, Cl moved substantially with Li diffusion, while Zr and O moved slightly during AIMD simulations by extracting MSD of each atom at 380 K (Supplementary Fig. 13b). These results highlight that the deficient Li1.75ZrCl4.75O0.5 may have Cl diffusion which may need further in-depth analyses.Fig. 5: Thermal cycle, MSD, Arrhenius plots and pair correlation function for amorphous structure.a Thermal cycle indicating heat-quench method for the formation of amorphous structure for Li1.75ZrCl4.75O0.5. The red part indicates the heating process to completely melt the structure at 1000 K, blue portion represents cooling with constant simulation time, and green line indicates the final equilibrated structure at 300 K. b MSD of Li extracted for amorphous Li1.75ZrCl4.75O0.5 at different temperatures, c Arrhenius plots of diffusivity at different temperatures (340, 360, and 380 K) to compute extrapolated Li ionic conductivity at room temperature for amorphous Li1.75ZrCl4.75O0.5. Pair correlation function for amorphous structure extracted from AIMD simulations for d Li–Li pair and e Li–Cl pair at each temperature.To check the effect of cell size and large number of atoms on the transport mechanism and Li ionic conductivity in an amorphous structure, a supercell of 2 × 1 × 2 with 96 atoms having lattice parameters of a = 22.48571 Å, b = 11.24285 Å, c = 13.01690 Å is used as a case study with same computational techniques with different temperature ranges and simulation time to melt the structure at 2000 K (for more detail see Supplementary note 3). The 96 atoms model exhibits fast Li ion diffusion at 300 K with Li-ion diffusion coefficients of 7.92 × 10−6 cm2/s and smoother than 48 atoms model whose diffusion coefficient is 7.92 ×10-6 cm2/s (Supplementary Fig. 14). The structural configuration at different temperatures for both models is presented in Supplementary Figs. (15−17) during the formation of amorphous phase. Through an extended simulation, it has been noted that the MSD of Cl is decreasing, suggesting that Cl is not actively involved in the diffusion process as a chloride ion. Rather, the motion of the Cl atom seems to be confined to rotation and oscillation around its periodic and mean position (Supplementary Fig. 18).To confirm the formation mechanism of the amorphous phase, PCF g(r) of Li1.75ZrCl4.75O0.5 was extracted at melting, cooling and 340 K temperature. The PCF of Li-Li pair at 1000 K (dark yellow line) indicates hyperbolic shape and no significant straight peak, which indicates that the structure is completed melted (Fig. 5d). After cooling, Li pair have different peaks indicating strong bonding length of Li with other atoms (for detail see Supplementary Fig. 19), while at 340 K Li pair appears with first peak at 3.97 Å and second high peak at radius of 6.15 Å, indicating fast diffusion and loosely bonding of Li. The inset Figure illustrates that the bonding distance between two Li atoms aligns with the PCF first peak, while the bonding length between Li and Cl falls within the range of 2.30−2.72 Å. Figure 5e shows that Li-Cl pair has the same bonding radius at each temperature (melting and cooling), with the first highest peak at the radius of ≈2.5 Å, that can be seen in the inset Figure where bonding length between Li and Cl varies in the range 2.39−2.74 Å, however the peak height decreases with temperature, showing Cl may have diffusion along Li.