Atomic-scale studies of garnet-type Mg3Fe2Si3O12: Defect chemistry, diffusion and dopant properties

Materials with low cost, environmentally benign, high structural stability and high Mg content are of considerable interest for the construction of rechargeable Mg-ion batteries. In the present study, atomistic simulations are used to provide insights into defect and diffusion properties of garnet type Mg3Fe2Si3O12. Calculations reveal that the Mg–Fe anti-site defect cluster (0.44 eV/defect) is the lowest energy intrinsic defect process. Three dimensional Mg-ion migration pathway with the activation energy of 2.19 eV suggests that Mg-ion diffusion in this material is slow. Favourable isovalent dopants are found to be Mn2þ, Ga3þ and Ge4þ on the Mg, Fe and Si sites respectively. While the formation of Mg interstitials required for the capacity is facilitated by Al doping on the Si site, Mg vacancies needed for the vacancy assisted Mg-ion diffusion are enhanced by Ge doping on the Fe site. The electronic structures of favourable dopants are calculated and discussed using density functional theory.


Introduction
The continued depletion of fossil fuel resources has led to the research activity to find alternative energy storage systems such as batteries and fuel cells. Rechargeable Li-ion battery is currently being considered as a potential energy source for many portable electronic devices [1][2][3][4][5]. As the global demand for large scale applications such as electrical vehicles increases rapidly, battery systems based on other monovalent cations such as Na þ [6,7] or divalent cations such as Mg 2þ [8], Ca 2þ [9] and Zn 2þ [10] are being considered as alternatives to Li-ion batteries. Batteries based on divalent cations are expected to provide high energy density for large-scale applications as they can undergo beyond one-electron redox reactions.
Mg-ion batteries are of considerable interest for the next generation of energy storage due to low cost, high energy density and high safety [11][12][13][14]. However, the intercalation of Mg 2þ ions in the electrode materials is expected to be difficult due to the strong ionic interaction between Mg 2þ ions and the host lattices. This problem is now being partly resolved by synthesizing electrode materials at the nano-scale [15] and reducing the shielding of Mg 2þ ions in the lattice [16]. A limited number of electrode materials including MgFeSiO 4 [17], MgMn 2 O 4 [18], Mg 6 MnO 8 [19], V 2 O 5 [20] and MnO 2 [21] have been so far examined.
Naturally occurring minerals have been considered for potential application as electrode materials in Mg and Ca-ion batteries [19,[22][23][24].
Recently MgFeSiO 4 which belongs to olivine-type mineral was examined as a promising cathode material for Mg-ion batteries [23]. Mg-based rock-forming silicate garnet Mg 3 Fe 2 Si 3 O 12 [25] is a candidate material for assessing its suitability as an electrode material for numerous reasons. This material consists of three Mg ions per formula unit and (SiO 4 ) 4units offering structural stability via strong Si-O bonds. Furthermore, silicon and iron are relatively safe, abundant and cheap. To the best of our knowledge, Mg 3 Fe 2 Si 3 O 12 was not studied either experimentally or theoretically to examine the defect, dopant, transport and electronic properties. As defect and transport properties of Mg 3 Fe 2 Si 3 O 12 are of greater importance in assessing its suitability in batteries, in this study, computer modelling based on the interatomic potentials is carried out using established lattice energy minimisation technique to provide information about the lattice properties, defect formation energies, Mg-ion diffusion pathways and solution of dopants. This method has been successfully applied to a wide range of energy materials [26][27][28][29]. Additionally, density functional theory (DFT) based energy minimisation techniques was applied to investigate the electronic structures of doped configurations.

Computational methods
Defect calculations were performed using a well-established classical simulation code GULP (General Utility Lattice Program) [30]. This code uses the Born model description of an ionic crystal lattice. Lattice energy is defined as the sum of long-range attractive forces (Coulomb) and short-range forces [repulsive (Pauli) and attractive (van der Waals)]. Buckingham potentials were used to model the short range interactions. Ionic positions and lattice parameters were relaxed using the Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm as implemented in the GULP code [31]. All defect calculations were performed with a supercell consisting of 704 atoms. Point defects were modelled using the Mott-Littleton method [32]. Two adjacent Mg vacancy sites were considered as initial and final configurations for the migration of Mg-ion. Activation energy was defined as the energy difference between the vacancy formation energy and the maximum energy along the diffusion path. As defects are at dilute limit with full charge, there will be an overestimation in the defect enthalpies. However, relative energies and trends will be consistent. From a thermodynamical point of view, the defect parameters (i.e. formation energies) can be defined via the comparison of the defective (real) crystal to an isobaric (or isochoric) non-defective (ideal) crystal. These defect formation parameters are interconnected through thermodynamic relations as outlined in previous studies [33,34]. In this ssudy, calculations correspond to the isobaric parameters for the migration and formation processes [35,36].
Electronic structure calculations were performed using a plane wave DFT code VASP (Vienna Ab initio Simulation Package) [37]. This code uses plane wave basis sets and PAW (Projected Augment Wave) pseudo-potentials [38]. A plane wave basis set with a cut-off of 500 eV was used in all calculations. Exchange-correlation effects were treated using generalized gradient approximation (GGA) as implemented by Perdew, Burke and Ernzerhof (PBE) [39]. The conjugate gradient (CG) algorithm was used to minimize the energy. A supercell containing 160 atoms was used. The force tolerance was less than 0.001 eV/Å in all relaxed configurations. A 2 Â 2 Â 2 Monkhorst Pack k-point mesh [40] which yielded 8 k-points was used. Short range dispersive attractive interactions were modelled using DFT þ D3 method implemented by Grimme et al. [41].

Results and discussion
3.1. Crystal structure of Mg 3 Fe 2 Si 3 O 12 Mg 2 Fe 2 Si 3 O 12 belongs to the garnet-type and crystallises in the cubic space group Ia3d (no. 230) with the lattice parameters [25]. The Mg 2þ ion forms 8-fold coordination with adjacent oxygens. Octahedral (FeO 6 ) and tetrahedral (SiO 4 ) units share their corners and form a three-dimensional network as shown in Fig. 1a. Energy minimisation calculations were performed to test the suitability of the classical pair potentials, pseudopotentials and basis sets used in this study. Table 1 reports the experimental and calculated lattice parameters, showing the  good agreement and supporting the validity of the parameters used in this study. The calculated maximum specific capacity of Mg 3 Fe 2 Si 3 O 12 is 349 mAhg -1 (see supplemenatry information for the calculation formula used) inferring the importance of this material for future application magnesium ion batteries. However, practical difficulties do not allow to provide maximum specific capacity and the real quantity should be calculated from the experiments.

Energetics of intrinsic defect processes
Classical simulations were applied to gain insight into the formation of point defects and key defect processes such as Schottky, Frenkel and anti-site. These defect processes are important as they can influence the electrochemical properties of Mg 3 Fe 2 Si 3 O 12 . The point defects (vacancies and interstitials) were combined to calculate Schottky and Frenkel defect formation energies. In order to write complete Schottky reaction equations (eqs (6)-(8)), appropriate lattice energies were calculated (refer to  Table S3 for lattice phases together with their lattice energy per formula). The following equations as written using Kr€ oger-Vink notation [42] describe those defect processes in detail.
Fe antisite ðisolatedÞ : In Fig. 2, we report the calculated reaction energies per defect for all processes. The lowest energy defect process is for the Mg-Fe anti-site defect cluster. Other anti-site defect cluster energies are also calculated to be low. This indicates a small amount of cation mixing can be found in this material. Many experimental and theoretical studies identified this defect at different experimental conditions and during cycling the batteries [43][44][45]. Schottky disorder is a process in which charge-balancing populations of vacancies are created and then displaced to the surface of the crystal (equation (5)). Among other defect processes, the MgO-Schottky is the lowest energy process with the defect energy of 3.28 eV/defect. This process will enhance the formation of V '' Mg and V O at high temperatures. The oxygen and magnesium Frenkel energies are calculated to be 3.44 eV/defect and 4.09 eV/defect respectively and they differ by 1 eV from the MgO-Schottky ensuring the possibility of MgO loss in this material at high temperatures. The Fe Frenkel, Si Frenkel and SiO 2 Schottky energies are highly endothermic suggesting that they will not be present at significant concentrations.

Diffusion of Mg-ions
In this section, we discuss the diffusion of Mg 2þ ions in order to assess the suitability of this material as a cathode in Mg-ion batteries. Performance of a battery is partly dependent on the materials with high ionic conductivity. Slow diffusion will degrade the battery performance during the intercalation. As diffusion depends on the crystal structure, examining diffusion pathways together with activation energies in Mg 3 Fe 2-Si 3 O 12 is necessary. Here we use classical simulation to calculate the vacancy assisted migration of Mg 2þ ions. In a previous study [46], this methodology has proven to be efficient. For example, there was an excellent agreement between the calculated Li-ion migration pathway and that observed in the neutron diffraction experiment in LiFePO 4 [46,47].
We identified a potential Mg-Mg hop with the jump distance of 3.54 Å. Importantly, this hop is connected to form a three-dimensional long-range migration path as shown in Fig. 1b. Activation energies calculated for linear and curved pathways are 2.78 eV and 2.19 eV respectively (refer to Fig. 1c). In the linear pathway, Mg 2þ ions migrate in a straight line and this pathway does not guarantee the lowest activation energy. We allowed the migration ion in different directions as described in our [26,28,29] and previous simulations [27,47] to attain the lowest

Solution of dopants
Substitution of suitable dopants is an effective way of stabilizing a material in order to prevent unfavourable transformations and enhance its electronic properties. The most favourable dopants can be of interest in the experimental investigations to tune the properties of Mg 3 Fe 2-Si 3 O 12 . A wide range of cations were substituted and their solution energies were calculated. Appropriate charge compensating defects and lattice energies were combined to construct the defect equations. The results are discussed in terms of solution energy, cation size and local environment of dopant atoms. Buckingham potentials used for dopants are given in the electronic supplementary information (ESI)(refer to Table S1). Formulae used to calculate defect energies and solution energies are given in the ESI (refer to Table S2). Table S3 in the ESI provides the information of lattice phases considered together with lattice energies per formula unit. Bond distances and bond angles calculated around the favourable dopants using DFT and classical simulations were in good agreement with a strong validation of the results.

Divalent dopants
First we considered the substitution of divalent cations (M ¼ Ni, Zn, Mn, Co, Ca, Sr and Ba) on the Mg site. The following reaction equation describes the doping mechanism in which there is no charge compensation. The calculated solution energies are shown in Fig. 3a.
The lowest solution enthalpy (-0.42 eV), and therefore highest solubility, is predicted for Mn 2þ . The negative solution enthalpy suggests that Mn 2þ would be a promising stabilizing agent. Notably, the solution enthalpy for Ca 2þ is also negative (À0.28 eV) meaning that it can also be tested experimentally. A very low endoergic solution enthalpy of 0.10 eV is calculated for Zn 2þ and Co 2þ suggesting that these two dopants are also candidates to be tested but would require a small amount of energy in the form of heat. The ionic radius of Mg 2þ in an eight coordination environment is 0.89 Å. Favourability of these four dopants can be due to their ionic radii (between 0.74 Å and 1.00 Å) which are close to the ionic radius of Mg 2þ ion. A small ionic radius difference (0.20 Å) between Ni 2þ and Mg 2þ results a slightly high endoergic (0.31 eV) solution entha;py. Solution enthalpy increases with ionic radius from Sr 2þ (0.77 eV) to Ba 2þ (2.70 eV). The highest solution enthalpy calculated for Ba 2þ suggests this dopant is very unfavourable on the Mg site.

Trivalent dopants
Thereafter, we considered a range of trivalent cations (M ¼ Al, Ga, In, Sc, Y, Gd and La) on the Fe and Si sites. Substitution on the Fe site does not require any point defects to charge compensate because of the þ3 charge on the Fe in Mg 3 Fe 2 Si 3 O 12 (refer to equation (16)). Solution enthalpies are plotted as a function of dopant size in Fig. 3b. The most favourable dopant is calculated to be Ga 3þ with the exothermic solution enthalpy of À0.03 eV. This can be due to the ionic radius of Fe 3þ (0.65 Å) closer to that of Ga 3þ (0.62 Å). The second most favourable dopant is Al 3þ and its solution enthalpy is 0.13 eV. Solution enthalpy increases with the ionic radius from Sc 3þ to La 3þ . The highest solution enthalpy is calculated for La 3þ (ionic radius 1.03 Å) meaning that La doping is unfavourable.
The Si site was also considered for doping to create additional Mg 2þ ions as charge compensating defects in this material (refer to equation (17) in what follows). This defect mechanism can be an efficient way of increasing the concentration of Mg in Mg 3 Fe 2 Si 3 O 12 in order to increase its capacity. In a previous experimental study, Co 3þ was doped on the Ru site in Li 2 RuO 3 and it was shown that doping enhanced the electrochemical reversibility of Li þ ions and the amount of Li þ in the lattice [54].
In Fig. 3d, we report the solution enthalpies for this defect mechanism. In all cases solution enthalpies are highly positive (>4.50 eV) showing the inefficacy of this process. However, the most favourable dopant Al 3þ might be worth trying experimentally. The lowest solution enthalpy for Al 3þ can be due to the small difference in their ionic radii. High endothermic process for this mechanism can be due to the high defect formation energy of quadruply charged Si.

Tetravalent dopants
Finally tetravalent dopants (M ¼ Ge, Ti, Sn, Zr and Ce) were considered for the substitution on the Si and Fe sites. Solution enthalpies calculated for the substitution on the Si site are reported in Fig. 3c. The following defect equation explains the mechanism involved in this process.
The results show that the most favourable solution energy (0.75 eV) is predicted for Ge 4þ. This is due to a small ionic radius difference between Ge 4þ (0.39 Å) and Si 4þ (0.26 Å). The solution enthalpy for Ti 4þ is calculated to be 4.29 eV meaning that it is an unfavourable dopant. The second favourable dopant is Sn 4þ with the solution enthalpy of 2.87 eV. Solution enthalpy gradually increases with increasing ionic radius. The highest enthalpy of solution for Ce 4þ (5.36 eV) suggests that doping would require external heat energy.
The formation of Mg vacancies required for the vacancy assisted Mgion migration was considered by doping tetravalent cations on the Fe site: Calculated energies of solution are reported in Fig. 3e. Germanium is identified as a promising dopant for this process as it has the lowest solution enthalpy of 0.79 eV. Solution enthalpy gradually increases with ionic radius and the highest positive enthalpy value is noted for Ce 4þ meaning that doping is unlikely to occur at normal temperatures.

Mg-ion diffusion in the presence of dopants
Here we calculate the activation energies for Mg-ion diffusion in the presence of Al 3þ on the Si site, Ge 4þ on the Si site and Ga 3þ on the Fe site. Energy profile diagrams for the local Mg hop are shown in the ESI (refer to Fig. S1). There is a reduction (by 0.07 eV) in the activation energy for the doping of Al 3þ on the Si site. In the case of doping of Ge 4þ on the Si site, activation energy is reduced by 0.09 eV. Doping of Ga 3þ on the Fe site results a reduction in the activation energy by 0.11 eV. In all cases, Mg-Mg distances have been slightly altered compared to the distances present in the un-doped crystal structure. This perturbation in the distances are due to the charge or ionic radius mismatch between Al 3þ and Si 4þ , Ge 4þ and Si 4þ and Ga 3þ and Fe 3þ . Current simulation shows that doping of Al 3þ and Si site would simultaneously increase the concentration of Mg 2þ ions in the lattice and slightly reduces the activation energy of the Mg-ion migration.

Electronic structure of doped Mg 3 Fe 2 Si 3 O 12
Here we use DFT simulations to examine the electronic structures of doped configurations and discuss the results calculated for the most favarouble dopants (Mn 2þ , Ga 3þ, Al 3þ and Ge 4þ ) as discussed in section 3.4.
The DOS plot calculated for bulk Mg 3 Fe 2 Si 3 O 12 is shown in Fig. 4a. It is shown that the ferromagnetic configuration of bulk Mg 3 Fe 2 Si 3 O 12 exhibits metallic behaviour. Doping Mn does not alter the electronic structure much as shown in the DOS plot (refer to Fig. 4b). Peaks associated with the 3d states of Mn appear near the Fermi level (refer to Fig. 4c) Fig. 4e.
In the case of Ga, the Fermi level is slightly altered but maintaining metallic behaviour (refer to Fig. 4f). There is no significant contribution to the DOS plot from Ga near the Fermi level. States associated with Ga appear in the conduction band (refer to Fig. 4g). Ga-O bond distances are longer only by~0.05 Å than the Fe-O bond distances (1.92 Å-1.94 Å). This is due to a very small diference in their ionic radii.
Doping Ge on the Si site has only a minor effect on the DOS plot (refer to Fig. 4j). The doped configuration is metallic. A very small peak that belongs to 3s and 3p states of Ge appear just below the Fermi level (refer to Fig. 4k). The Si-O bond distance in the undoped Mg 3 Fe 2 Si 3 O 12 is 1.65 Å. The Ge-O bond distance is longer than (by 0.15 Å) than the Si-O bond distance. This is due to the larger ionic radius of Ge 4þ (0.39 Å) than that of Si 4þ (0.26 Å).
Doping Al on the Si site shifts the Fermi energy level only by 0.07 eV (refer to Fig. 5b). The doped configuration is metallic. Additional peaks arising from 3s and 3p states of Al appear near the Fermi level (refer to Fig. 5c). Calculated Al-O bond distance (1.77 Å) is longer than the Si-O bond distance (1.65 Å). This is due to the larger ionic radius and smaller charge density of Al 3þ than that of Si 4þ .
There is no shift in the Fermi energy upon doping of Ge on the Fe site (refer to Fig. 5f). Very small peaks associated with Ge appear near the Fermi level (refer to Fig. 5g). Metallic behavior is still maintained. A small difference in the Ge-O bond distance (1.91 Å) is noted comprared to Fe-O bond distances. This is due to the small difference in the ionic radius and charge mismatch. Constant charge density plot of Ge-doped on the Fe site in Mg 3 Fe 2 Si 3 O 12 is shown in Fig. 5h.

Conclusions
In the present study, we used classical and DFT simulations to examine the defect, diffusion, dopant and electronic properties of garnet type Mg 3 Fe 2 Si 3 O 12 . The Mg-Fe anti-site defect cluster is the lowest energy intrinsic defect process. Mg-ion diffusion in this material is three-dimensional with a high activation energy of 2.19 eV and therefore it is low. The Mg, Fe and Si sites prefer to be doped by isovalent Mn 2þ , Ga 3þ and Ge 4þ ions respectively. Capacity of this material can be increased by the formation of additionnal Mg in the form of interstitials by doping of Al 3þ on the Si site. The formation of Mg vacancies required for the vacancy assisited migration can be facilitated by the doping of Ge 4þ on the Fe site. Doping reduced the diffusion barrier only by 0.1 eV suggesting that it would be difficult to make this material suited ion conductor.

Funding
This research was financially supported by European Union's H2020 Programme under Grant Agreement no 824072-HARVESTORE.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.