Computational investigation of the Mg-ion conductivity and phase stability of MgZr₄(PO₄)₆

Abstract

Solid electrolyte materials exhibiting high Mg-ion conductivity are required to develop Mg-ion batteries. In this study, we focused on a Mg-ion-conducting solid phosphate based electrolyte, MgZr₄(PO₄)₆ (MZP), and evaluated the ionic conductivity of NASICON-type and β-iron sulfate-type MgZr₄(PO₄)₆ structures via density functional theory calculations. The calculations suggest that the migration energy of Mg is 0.63 eV for the NASICON-type structure and 0.71 eV for the β-iron sulfate-type one, and the NASICON-type structure has higher ion conductivity. Although the NASICON-type MZP structure has not been experimentally realised, there is only an energy difference of 14 meV per atom with respect to that of the β-iron sulfate-type structure. Therefore, in order to develop a synthesis method for the NASICON-type structure, we investigated pressure- and temperature-dependent variations in the free energy of formation using density functional perturbation theory calculations. The results suggest that the formation of the NASICON-type structure is disfavoured under the 0–2000 K and 0–20 GPa conditions.

---

Introduction

Lithium ion batteries (LIBs) are widely used in portable electronic devices, such as mobile phones and laptop computers, owing to their high energy density and long cycle life. In recent years, the use of (hybrid) electric vehicles (EVs) has gradually increased to address the environmental and energy issues.¹ Therefore, LIBs are also attracting attention as secondary batteries for on-vehicle use. However, the following two major technical drawbacks of current LIBs need to be overcome: (i) insufficient driving range of the vehicle compared to the conventional petro-powered vehicle due to the small energy density and (ii) safety concerns arising from inflammable organic electrolytes. In this respect, developing all solid-state Mg-ion batteries (MIBs) is believed to be one of the ultimate solutions to resolve the above two issues.² The Mg ion is divalent and can carry two electrons, so that the energy density is expected to double. All solid-state battery technology also contributes to the fabrication of a bipolar-type battery, thus enabling increased energy density, and replacement of organic electrolyte by an inflammable solid electrolyte can solve the safety concerns.

However, the generally poor Mg diffusivity in the solid state prevents the realisation of MIB-based devices, since divalent Mg ions strongly interact with counter anions. Therefore, it is necessary to understand the factors affecting the Mg-ion conductivity and find fast Mg-ion conductors to develop future all solid-state MIBs.

As solid electrolyte materials, oxide-based Na super ionic conductor (NASICON)-type structures are attractive for LIBs because NASICON-type ion conductors show high ion conductivity in both Na^3–9 and Li^10–17 systems. So far, one of the best Mg-ion conductors is rhombohedral NASICON-type (Mg_0.1Hf_0.9)_4/3.8Nb(PO₄)₃ (Fig. 1(a and b)), in which the ionic conductivity of Mg²⁺ is 2.1 × 10⁻⁶ S cm⁻¹ at 573 K and its activation energy is 0.68 eV.¹⁸ This compound consists of (Nb,Hf)O₆ octahedra and PO₄ tetrahedra alternately bridged at the two-coordinated edge O atoms, and Mg ions migrate through the interstitial sites of the corner-shared network. This local bonding rule leads to another type of structure with a different network topology, known as the monoclinic β-iron sulfate-type structure, and MgZr₄(PO₄)₆ belonging to this structure type (Fig. 1(c and d)) also exhibits relatively high Mg-ion conductivity.^19,20 An example of a material with this type of structure is Mg_0.7(Zr_0.85Nb_0.15)₄(PO₄)₆. The ionic conductivity of this compound is 1.1 × 10⁻⁷ S cm⁻¹ at 573 K and the activation energy is 0.92 eV,²¹ and Mg-ion conduction of this compound at 350 °C has been recently confirmed using non-blocking Mg metal electrodes.²² Interestingly, both the NASICON-type and β-iron sulfate-type structural frameworks have a polymorphic relationship represented by A_xB₄(XO₄)₆ frameworks. In case of Na ion migration, the NASICON-type structure has been reported to show higher conductivity than that of the β-iron sulfate-type one for the composition of Na₃Zr₂Si₂PO₁₂.⁴ The lower activation energy for the Mg ion conduction in the NASICON-type (Mg_0.1Hf_0.9)_4/3.8Nb(PO₄)₃ structure than in the β-iron sulfate-type Mg_0.7(Zr_0.85Nb_0.15)₄P₆O₂₄ structure has been attributed to the higher structural symmetry of the former structure.¹⁵ NASICON-type MgZr₄(PO₄)₆ has not been synthesised primarily because the stability of the NASICON-type structure is lower for small interstitial cations like Mg²⁺ ions.¹⁹

Fig. 1 (a) and (b) show the NASICON-type structure of MZP and (c) and (d) show the β-iron sulfate-type structure of MZP with the ZrO₆ octahedra depicted in green, PO₄ tetrahedra depicted in lavender, O atoms in red, and Mg atoms in orange.

However, the relationship between the structure and ionic conductivity is not yet clear for the case of Mg ions, because there are no experimental ionic conductivity data of both structures with the same composition. In this study, we evaluate the ionic conductivity for both the NASICON-type and β-iron sulfate-type structures with the stoichiometric composition of MgZr₄(PO₄)₆ (MZP) using first-principles molecular dynamics (FPMD) study based on density functional theory (DFT). So far, MZP has been experimentally prepared only in the β-iron sulfate-type structure, and no report on the synthesis of the NASICON-type structure is known. We also investigated the phase stability of NASICON-type and β-iron sulfate-type structures under various temperature and pressure conditions.

Experimental section

Computational details

Structure relaxation parameters including the lattice parameters and internal coordinates were calculated by first-principles DFT calculations. In detail, a Vienna ab initio simulation package (VASP)^23–26 based on DFT²⁷ with a projector augmented-wave (PAW)^28,29 method and a plane-wave basis set was used. Further, we used a generalized gradient approximation (GGA)-type exchange–correlation functional developed by Perdew, Burke, and Ernzernhof and later modified for solid materials (PBEsol).³⁰ The cutoff energy for the plane-wave basis was set to 500 eV and k-point resolution was set to ∼1000 (i.e., N_x, N_y, and N_z, the number of grids in k_x-, k_y-, and k_z-directions of the reciprocal space were set to satisfy N_x × N_y × N_z × N_atom ≈ 1000, where N_atom is the number of atoms in each unit cell).

The NASICON-type structure was constructed by replacing Li ions in the structure of LiZr₂(PO₄)₃ taken from the inorganic crystal structure database (ICSD) (ID: 92250, original data reported by Catti et al.)³¹ with Mg ions and vacancies, because there are no reports of the experimental synthesis. The crystal structure input for β-iron sulfate-type MZP was extracted from the ICSD (ID: 250452, original data reported by Gobechiya et al.).³²

Table 1 lists DFT-derived relaxed lattice parameters, cell volume, and calculated total electron energies for (a) NASICON-type and (b) β-iron sulfate-type MZP. In the NASICON-type structure, three Mg ions and vacancies reside at the 6b and/or 36f Wyckoff position (36f sites are split around the 6b site). The most stable Mg/vacancy arrangement was determined using a genetic algorithm (GA), as described in our previous paper.^33–35 Mg-ion occupation at 6b sites is the most stable after GA optimisation, and the energy difference among the GA-derived structures (2400 samples) is within 60 meV per atom.

Table 1 DFT-derived relaxed lattice parameters, cell volume, and calculated total electron energies for the most stable structure of (a) NASICON-type and (b) β-iron sulfate-type MZP

	a/Å	b/Å	c/Å	α/°	β/°	γ/°	Cell volume per atom/Å^3 per atom	Electron energy per atom/eV per atom
(a)	8.91	8.98	22.4	89.80	89.99	119.99	14.8	−8.041
(b)	12.52	8.97	8.95	90.00	90.49	90.00	14.4	−8.055

---

In contrast, Mg ions and vacancies occupied the 4e sites of the Wyckoff position at 1 : 1 molar ratio in the β-iron sulfate-type MZP. Three configurations of Mg/vacancy arrangements were calculated for the unit cell of Mg₂Zr₈P₁₂O₄₈, and the most stable one was chosen for following calculations. The volumetric and energetic differences among the three configurations are 0.19 ų per atom and 12 meV per atom, respectively, indicating a small dependence on the Mg arrangement.

The obtained Mg/vacancy configurations for NASICON-type and β-iron sulfate-type MZP structures and their lattice parameters were thereafter used as structural inputs for following calculations, unless specified otherwise.

As listed in Table 1, the β-iron sulfate-type MZP is slightly more stable than the NASICON-type one by only 14 meV per atom, and the cell volume per atom difference between β-iron sulfate-type MZP (14.4 ų per atom) and NASICON-type (14.8 ų per atom) is small. Therefore, both structures are reasonably stable and are expected to be good candidates for comparing the Mg-ion diffusivity.

FPMD calculations were performed to compare the Mg-ion conductivity of the NASICON-type and the β-iron sulfate-type MZP structures. Initially, all the Mg atoms were located in the 6b sites for the NASICON-type structure and 4e sites for the β-iron sulfate-type structure.

The supercell structure including eight MZP units (i.e., Mg₈Zr₃₂P₄₈O₁₉₂) was prepared to evaluate the Mg-ion conductivity via FPMD simulations. The lattice constants of the supercell of the NASICON-type and β-iron sulfate-type MZP are a = 12.54 Å, b = 17.62 Å, c = 17.64 Å, α = 89.99°, β = 90.01°, γ = 89.00° and a = 12.52 Å, b = 17.94 Å, c = 17.90 Å, α = 90.00°, β = 90.49°, and γ = 90.00°, respectively. The lattice vectors a_sc, b_sc, and c_sc in the supercell are described by a linear combination of lattice vectors a_p, b_p, and c_p in a primitive cell of the NASICON-type or β-iron sulfate-type structure, as shown below:

(a_scb_scc_sc) = (a_pb_pc_p)M_sc

where, the supercell matrix M_sc is expressed as follows:

FPMD simulations were performed for NASICON-type and β-iron sulfate-type structure in the temperature range of 1573 to 1973 K with increments of 50 K for 100 ps. The DFT calculation settings used were equivalent to those of the structure relaxation calculations except for the cutoff energy and k-point grids as follows: the cutoff energy for the plane-wave basis was set to 350 eV and 1 × 1 × 1 k-point grid (only Γ point) was employed to reduce the computational cost. For MD simulations, the time step was set to 1 fs, and the NVT canonical ensemble was adopted using a Nosé thermostat.³⁶ The mean square displacements (MSD) of all elements were calculated, and then the diffusion coefficients of Mg-ion conduction and corresponding migration energies were determined. Above FPMD studies have been performed in our previous studies, showing reasonable accordance between the experimental and the computational results.^10,11,37

We evaluated the thermodynamic stability of the MZP polymorphs by first principles DFT calculations. Temperature- and/or pressure-dependent phase stability was also evaluated using DFT-based phonon calculations or cell-volume constraint calculations, respectively. We used the same lattice as that used for FPMD simulations. Structural relaxation calculations were performed for a structure in which the lattice constant was varied from 95 to 103% under the condition of constant lattice volume. The energy–volume relationship is fitted with the Murnaghan equation of state:^38,39

where, B₀ is the bulk modulus at zero pressure, B′ is its first derivative, E₀ is the minimum energy, and V₀ is the volume at the minimum energy. In addition, the bulk modulus B can be expressed as,

B = −Vdp/dV (2)

where, V and p are the volume and pressure, respectively. The enthalpy–pressure diagram was constructed by combining eqn (1) and (2).

Phonon frequencies and their contribution to free energy changes were calculated using the PHONOPY code.⁴⁰ Density functional perturbation theory (DFPT),^41,42 also implemented in VASP, was used to estimate the real-space force constants of the supercells. The other calculation conditions were the same as those used for the structure relaxation calculations.

From these calculations, the free energy was calculated by changing the temperature from 0 to 2000 K and the pressure was changed from 0 to 20 GPa.

Results and discussion

Fig. 2 shows the evolution of the time-averaged mean square displacement (MSD) for constituent ions in MZP with β-iron sulfate-type structure. The MSD of Mg ions increases linearly with the MD step (time), indicating that Mg²⁺ ions diffuse over the lattice. In contrast, the MSD profiles of the other ions remains constant (<1.1 Ų). Therefore, Zr, P, and O ions remain near their original sites during the MD run, and these MSD values represent the magnitude of thermal vibrations. Similar MSD profiles were obtained for the NASICON-type structure, in which only Mg ions are found to diffuse in the lattice.

Fig. 2 MSD plots of Mg, Zr, P, and O trajectories in β-iron sulfate-type MZP at 1973 K.

Diffusion path of Mg ions at 1573 K and 1973 K, i.e. both the lowest and the highest simulation temperatures, are visualised through isosurfaces of the Mg-ion probability density distribution from the FPMD simulation for (Fig. 3(a and b) and ESI Fig. S1(a)†) the NASICON-type and (Fig. 3(c and d) and ESI Fig. S1(b)†) β-iron sulfate-type structures. We confirmed the formation of a 3-dimensional pathway for both structures. Moreover, metastable Mg-ion sites are indicated in both (Fig. 3b) NASICON-type and (Fig. 3d) β-iron sulfate-type structures by varying the isosurface level (Li ion probability density). In Fig. 3b and d, the stable and metastable sites are marked with orange and blue circles, respectively. On one hand, the metastable Mg sites in β-iron sulfate-type structure are linked to three neighbouring 4e stable Mg/vacancy sites with a straight diffusion path. On the other hand, the metastable Mg sites (18e sites) are connected to two neighbouring 6b stable sites, which is consistent with the reported Li migration pathway in LiZr₂(PO₄)₃ and its derivative.⁴³

Fig. 3 Trace of Mg atom in FPMD simulation at 1973 K. (a), (b) and (c), (d) show the NASICON-type and β-iron sulfate-type structures, respectively. In (b) and (d), the stable and metastable sites of Mg are denoted in orange and blue, respectively.

Diffusion coefficients were extracted through a simple linear fitting with the following expression:

〈x²〉 = 2nDt (3)

where, n, D, and t represent the dimension of Mg-ion diffusion, ionic diffusion coefficient, and time, respectively. The ionic conductivity of Mg is estimated by the Nernst–Einstein relationship:where, R, T, z, F, and c represent the gas constant, temperature, ionic valence, Faraday constant, and ionic concentration, respectively. Since number of Mg²⁺ hopping events are limited within the simulation time (100 ps) at ambient temperature according to Boltzmann statistics and attempt frequency, high temperature simulation (>1573 K) is necessary to evaluate migration energy in sufficient accuracy. Therefore, low temperature conductivities are extrapolated from high temperature data as reported in our previous FPMD calculations.^10,11,37

The Arrhenius plots for the NASICON-type and the β-iron sulfate-type structures is shown in Fig. 4, and the migration energies of Mg were evaluated by a simple linear fitting scheme. Higher Mg-ion conductivity and lower migration energy are derived for the NASICON-type structure compared to those of the β-iron sulfate-type structure.

Fig. 4 Arrhenius plots of the temperature-dependent Mg-ion diffusion coefficients fitted with a straight line.

The migration energy of Mg in the β-iron sulfate-type structure obtained in this study (0.71 eV) is different from the experimental value (1.6 eV),²¹ which may owe to inclusion of grain boundary conduction for experimental results. The β-iron sulfate-type structure has a higher migration energy than the NASICON-type structure (0.63 eV) in present computations. This tendency agrees that the migration energy of β-iron sulfate-type Mg_0.7(Zr_0.85Nb_0.15)₄(PO₄)₆ (0.92 eV)²¹ is higher than that of NASICON-type (Mg_0.1Hf_0.9)_4/3.8Nb(PO₄)₃ (0.68 eV).¹⁸ This result suggests that the NASICON-type structure has better ion conductivity than the β-iron sulfate-type one.

Since NASICON-type MZP has higher Mg-ion conductivity, the pressure- and temperature-dependence of the phase stability was investigated for NASICON-type and β-iron sulfate-type structures. Fig. 5 displays the difference between the phase stability of β-iron sulfate-type structure and NASICON-type one, ΔE, as functions of temperature and pressure. ΔE is described as follows:

ΔE = E(β-iron sulfate-type) − E(NASICON-type) (5)

where, E(β-iron sulfate-type) and E(NASICON-type) represent the free energies of the corresponding structures at a given pressure and temperature. Since ΔE is negative for all pressure and temperature conditions, the formation of the NASICON-type structure is thermodynamically disfavoured. However, the energy difference is only 14 meV per atom between two structures (at 0 K and at 0 atm), so that formation of NASICON-type structure may be feasible such as by doping techniques. Also, since the NASICON-type structure is stable under negative pressure condition, thin film formation technique using selection of appropriate substrate may stabilizes NASICON-type structure. In addition, we computed phase stability of MgTi₄(PO₄)₆ composition by DFT. NASICON-type structure is more stable than β-iron sulfate-type structure by 7 meV per atom, showing accordance with experimental formation of NASICON-type structure.⁴⁴ These results validate the present computational results.

Fig. 5 Temperature- and pressure-dependent free energies of β-iron sulfate-type and NASICON-type MZP.

Conclusions

First-principles DFT calculations were performed to clarify the difference between the Mg-ion conductivity in β-iron sulfate-type and NASICON-type MZP. The results of the FPMD calculations suggest that the NASICON-type structure has higher Mg-ion conductivity than the β-iron sulfate-type one. However, evaluation of the pressure- and temperature-dependent free energies by DFT calculations indicated that the formation of NASICON-type MZP is disfavoured. Therefore, the composition should be appropriately optimised to stabilize the NASICON-type structure.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

This work was partially supported by the ALCA-SPRING program from the Japan Science and Technology Agency (JST), a Grant-in-Aid for Scientific Research (18K19129) from the Japanese Ministry of Education Culture, Sports, Science and Technology (MEXT), and the MEXT program “Elements Strategy Initiative to Form Core Research Center” (since 2012). M. N. and Y. N. thank the support by the “Materials research by Information Integration” Initiative (MI²I) project of the “Support Program for Starting Up Innovation Hub” from the JST. Figures of crystal structures in this paper are drawn with VESTA.⁴⁵

References

1. R. Van Noorden, Nature, 2014, 507, 26–28.
2. H. Srour, L. Chancelier, E. Bolimowska, T. Gutel, S. Mailley, H. Rouault and C. C. Santini, J. Appl. Electrochem., 2016, 46, 149–155.
3. Y. Deng, C. Eames, L. H. B. Nguyen, O. Pecher, K. J. Griffith, M. Courty, B. Fleutot, J.-N. Chotard, C. P. Grey, M. S. Islam and C. Masquelier, Chem. Mater., 2018, 30, 2618–2630.
4. J. B. Goodenough, H.-P. Hong and J. A. Kafalas, Mater. Res. Bull., 1976, 11, 203–220.
5. H.-P. Hong, Mater. Res. Bull., 1976, 11, 173–182.
6. W. Zhou, Y. Li, S. Xin and J. B. Goodenough, ACS Cent. Sci., 2017, 3, 52–57.
7. F. Lalère, J. B. Leriche, M. Courty, S. Boulineau, V. Viallet, C. Masquelier and V. Seznec, J. Power Sources, 2014, 247, 975–980.
8. Y. Noguchi, E. Kobayashi, L. S. Plashnitsa, S. Okada and J. Yamaki, Electrochim. Acta, 2013, 101, 59–65.
9. A. G. Jolley, G. Cohn, G. T. Hitz and E. D. Wachsman, Ionics, 2015, 21, 3031–3038.
10. Y. Noda, K. Nakano, H. Takeda, M. Kotobuki, L. Lu and M. Nakayama, Chem. Mater., 2017, 29, 8983–8991.
11. Y. Noda, K. Nakano, M. Otake, R. Kobayashi, M. Kotobuki, L. Lu and M. Nakayama, APL Mater., 2018, 6, 060702.
12. X. Zhao, Z. Zhang, X. Zhang, B. Tang, Z. Xie and Z. Zhou, J. Mater. Chem. A, 2018, 6, 2625–2631.
13. P. Knauth, Solid State Ionics, 2009, 180, 911–916.
14. V. Thangadurai and W. Weppner, Ionics, 2006, 12, 81–92.
15. J. C. Bachman, S. Muy, A. Grimaud, H.-H. Chang, N. Pour, S. F. Lux, O. Paschos, F. Maglia, S. Lupart, P. Lamp, L. Giordano and Y. Shao-Horn, Chem. Rev., 2016, 116, 140–162.
16. X. Yao, B. Huang, J. Yin, G. Peng, Z. Huang, C. Gao, D. Liu and X. Xu, Chin. Phys. B, 2016, 25, 018802.
17. C. Cao, Z.-B. Li, X.-L. Wang, X.-B. Zhao and W.-Q. Han, Front. Energy Res., 2014, 2, 25.
18. S. Tamura, M. Yamane, Y. Hoshino and N. Imanaka, J. Solid State Chem., 2016, 235, 7–11.
19. K. Nomura, S. Ikeda, K. Ito and H. Einaga, Bull. Chem. Soc. Jpn., 1992, 65, 3221–3227.
20. S. Ikeda, M. Takahashi, J. Ishikawa and K. Ito, Solid State Ionics, 1987, 23, 125–129.
21. N. Imanaka, Y. Okazaki and G. Adachi, Electrochem. Solid-State Lett., 1999, 3, 327.
22. K. Kajihara, H. Nagano, T. Tsujita, H. Munakata and K. Kanamura, J. Electrochem. Soc., 2017, 164, A2183–A2185.
23. G. Kresse and J. Hafner, Phys. Rev. B: Condens. Matter Mater. Phys., 1993, 47, 558–561.
24. G. Kresse and J. Furthmüller, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 54, 11169–11186.
25. G. Kresse and J. Furthmüller, Comput. Mater. Sci., 1996, 6, 15–50.
26. G. Kresse and J. Hafner, Phys. Rev. B: Condens. Matter Mater. Phys., 1994, 49, 14251–14269.
27. P. Hohenberg and W. Kohn, Phys. Rev., 1964, 136, B864–B871.
28. G. Kresse and D. Joubert, Phys. Rev. B: Condens. Matter Mater. Phys., 1999, 59, 1758–1775.
29. P. E. Blöchl, Phys. Rev. B: Condens. Matter Mater. Phys., 1994, 50, 17953–17979.
30. J. P. Perdew, A. Ruzsinszky, G. I. Csonka, O. A. Vydrov, G. E. Scuseria, L. A. Constantin, X. Zhou and K. Burke, Phys. Rev. Lett., 2008, 100, 136406.
31. M. Catti and S. Stramare, Solid State Ionics, 2000, 136–137, 489–494.
32. E. R. Gobechiya, M. V. Sukhanov, V. I. Pet'kov and Y. K. Kabalov, Crystallogr. Rep., 2008, 53, 53–59.
33. G. Trimarchi and A. Zunger, Phys. Rev. B: Condens. Matter Mater. Phys., 2007, 75, 104113.
34. N. L. Abraham and M. I. J. Probert, Phys. Rev. B: Condens. Matter Mater. Phys., 2006, 73, 224104.
35. S. M. Woodley, P. D. Battle, J. D. Gale and C. R. A. Catlow, Phys. Chem. Chem. Phys., 1999, 1, 2535–2542.
36. S. Nosé, J. Chem. Phys., 1984, 81, 511–519.
37. R. Jalem, Y. Yamamoto, H. Shiiba, M. Nakayama, H. Munakata, T. Kasuga and K. Kanamura, Chem. Mater., 2013, 25, 425–430.
38. F. Birch, Phys. Rev., 1947, 71, 809–824.
39. F. Birch, J. Appl. Phys., 1938, 9, 279–288.
40. Togo, Phonopy, https://atztogo.github.io/phonopy/, (accessed 5 December 2018).
41. K. Refson, P. R. Tulip and S. J. Clark, Phys. Rev. B: Condens. Matter Mater. Phys., 2006, 73, 155114.
42. S. Baroni, S. de Gironcoli, A. Dal Corso and P. Giannozzi, Rev. Mod. Phys., 2001, 73, 515–562.
43. P. Padma kumar and S. Yashonath, J. Phys. Chem. B, 2001, 105(29), 6785–6791.
44. H. Takahashi and H. Takamura, Mater. Trans., 2012, 53, 932–93545.
45. K. Momma and F. Izumi, J. Appl. Crystallogr., 2011, 44, 1272–1276.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: 10.1039/c9ra00513g

---