Defect, transport, and dopant properties of andradite garnet Ca₃Fe₂Si₃O₁₂

Although lithium-ion batteries have dominated for years,^1–4 calcium-ion batteries (CIBs) are gaining attraction for large-scale storage applications due to cost, abundance, and enhanced safety.^5,6 Compared to other multivalent battery systems such as Mg-ion batteries,^7,8 the research activity of CIBs is in the early stage. CIBs are also expected to provide high energy and power density due to their two-electron reaction during intercalation. Furthermore, the electrode potential of Ca/Ca²⁺ is lower by 0.5 eV than that of Mg/Mg²⁺ and higher only by 0.1 eV than that of Li/Li⁺, implying higher cell voltage.⁹ 

Identification of suitable electrode materials is important in the development of CIBs. There are only a few published studies on electrode materials including layered compounds (e.g., V₂O₅),¹⁰ chevrel phases [e.g., CaMo₆X₈ (X = S, Se, and Te)],¹¹ and Prussian blue analogs [e.g., NaMnFe(CN)₆]⁹ for CIBs. Arroyo-de Dompablo et al.¹² recently carried out a joint experimental and theoretical study on CaMn₂O₄ polymorphs and concluded that full Ca extraction from these polymorphs introduces a substantial volume change with an average voltage of 3.1 V and a high Ca-ion migration barrier. In a recent theoretical study, Torres et al.¹³ considered a few Ca based minerals including Ca₃Mn₂(SiO₄)₃, CaMn(SiO₃)₂, and CaMn(CO₃)₂ to calculate their diffusion barriers, average voltages, and volume changes over Ca de-insertion.

Materials consisting of iron with structures formed by linking through polyhedral units such as PO₄, SiO₄, and SO₄, to provide structural stability, are of considerable interest in the development of electrode materials.^14–16 Andradite (Ca₃Fe₂Si₃O₁₂)¹⁷ is a rock-forming silicate garnet and consists of three Ca²⁺ ions per formula unit, leading to high theoretical capacity, and SiO₄⁴⁻ units offering structural stability via strong Si–O bonds. Furthermore, the low cost and availability of iron and silicon coupled with their non-toxicity also make this material ideal.

Computer modeling techniques based on pair-wise potentials are powerful tools to predict and understand the defects, ion transport, and dopant properties in solid-state oxide materials. A variety of oxide materials including battery materials have been examined using this technique.^18–22 In the present study, using inter-atomic pair potentials, we characterize the defect, diffusion, and dopant properties of Ca₃Fe₂Si₃O₁₂ on an atomistic level. Furthermore, density functional theory (DFT) allowed us to study the electronic structures of doped configurations.

Defect, diffusion, and dopant calculations were performed using the classical pair potential method as implemented in the GULP code,²³ which describes ionic interactions in the form of long-range Coulomb attractions and short-range attractions (electron–electron repulsion and attractive dispersion). Buckingham potentials (refer to Table S1) were used to describe short-range interactions. Cell parameters and ionic positions were relaxed using the Broyden–Fletcher–Goldfarb–Shanno (BFGS) algorithm.²⁴ Defects were modeled using the Mott–Littleton method.²⁵ Ca-ion migration was calculated by considering two adjacent Ca vacancy sites as initial and final configurations. Seven interstitial Ca ions were considered in a direct linear route, and they were allowed to relax. The difference between the vacancy formation energy and the maximum energy along the diffusion path is defined as the activation energy of migration. Relative energies and trends will be consistent although the present model assumes that ions are fully charged and defect concentrations are at a dilute limit. Thermodynamically, the defect parameters (i.e., migration energies) can be defined through the comparison of the defective crystal to an isobaric (or isochoric) non-defective crystal. These defect formation parameters are interlinked by thermodynamic relations.^26,27 The present calculations correspond to the isobaric parameters for the formation and the migration processes.^28,29

Electronic structure calculations are based on DFT using plane wave basis sets and projected augmented wave (PAW) potentials. A plane wave DFT code VASP (Vienna Ab initio Simulation Package)³⁰ which diagonalizes the Kohn–Sham equations was used. The calculations were performed using the generalized gradient approximation (GGA) as proposed by Perdew, Burke, and Ernzerhof (PBE).³¹ In all cases, a plane wave basis set with a cut-off of 500 eV and a 2 × 2 × 2 Monkhorst and Pack³² k-point mesh were used. Energy minimization calculations were performed using the conjugate gradient (CG) algorithm³³ and a force tolerance value of 0.001 eV/Å. Short-range attractive interactions (dispersion) were modeled using a DFT + D3 method as implemented by Grimme et al.³⁴ 

The crystal structure of Ca₃Fe₂Si₃O₁₂ at room temperature is cubic (space group $Ia3d¯$⁠, no. 230).¹⁷ Experimental lattice parameters reported by Novak and Gibbs.¹⁷ are a = b = c = 12.058 Å, α = β = γ = 90°, and V = 1753.18 ų. An eight-fold coordination is observed for Ca²⁺ cations. The Fe²⁺ ions reside in the FeO₆ octahedra. Tetrahedral (SiO₄) units share their corners with FeO₆ units in a three-dimensional network, as shown in Fig. 1. Energy minimization calculation was performed to relax both ionic positions and cell dimensions in order to obtain equilibrium lattice parameters. The results allowed us to check the quality of the Buckingham potentials used in the classical simulation and pseudopotentials together with basis sets used in the DFT simulation. There is good agreement between calculated and experimental values (refer to Table I).

Diffusion of ions in an ionic material is typically influenced by Schottky and Frenkel disorders, which can be calculated by combining point defects (vacancies and interstitials). These point defect energies can be difficult to determine via experimentation. The current classical simulation technique allowed us to calculate a series of point defects. The lowest energy point defects were then combined to calculate Frenkel and Schottky defect processes. Anti-site defects in which cations exchange their positions were also calculated as these defects have been observed in many oxide materials both experimentally and theoretically.^35–40 Here, we write reaction equations for all possible Schottky, Frenkel, and anti-site defects in Ca₃Fe₂Si₃O₁₂ using the Kröger–Vink notation,⁴¹ 

Figure 2 shows the resulting defect energies. The most favorable intrinsic disorder type is the Ca–Fe anti-site defect cluster (0.83 eV/defect). In this defect process, it is anticipated that there would be some Ca²⁺ ions on the Fe site and Fe³⁺ ions on the Ca site. In the isolated form of the anti-site defect, defects were modeled individually, and their energies were combined. Both defects were considered simultaneously in the cluster form. The energy difference between the isolated and cluster forms is considered as the binding energy (−0.19 eV). Notably, other anti-site defect cluster energies (Ca–Si and Fe–Si) are lower than the Schottky and Frenkel defect energies. Among other defect processes, the CaO Schottky is the lowest energy process, with a defect energy of 3.73 eV/defect. This process will facilitate the formation of $VCa″$ and $VO⋅⋅$ at high temperatures. The oxygen Frenkel energy is calculated to be 3.98 eV/defect, with a deviation of only 0.25 eV from the CaO Schottky. However, the Ca Frenkel energy is calculated to be ∼5.00 eV, inferring that the formation of Ca vacancies is not significant. Other Frenkel and Schottky energies are highly endoergic, suggesting that they will not be present at significant concentrations.

The diffusion of Ca²⁺ ions is discussed in this section. A material with high ionic conductivity generally exhibits a high rate performance and is suitable as an electrode for use in devices such as batteries. Classical simulation allowed us to calculate the Ca-ion diffusion pathways and the corresponding activation energies in Ca₃Fe₂Si₃O₁₂. In general, diffusion pathways are difficult to examine in the experiment, particularly for complex oxides such as Ca₃Fe₂Si₃O₁₂. The current classical simulation technique has been successfully applied in previous studies to calculate ion migration pathways together with activation energies.^42,43 For example, the Li-ion migration pathway calculated in LiFePO₄ by Fisher et al.⁴⁴ was in excellent agreement with the neutron diffraction experiment.⁴⁵ 

A potential Ca–Ca hop with a jump distance of 3.70 Å was identified. Many of these identical hops were then connected to form a three-dimensional long-range migration path, as shown in Fig. 3. The activation energy for this hop is calculated to be 2.63 eV, and Ca-ions migrate in a curved pathway. High activation energy of Ca²⁺ ion migration in Ca₃Fe₂Si₃O₁₂ shows that the Ca-ion diffusion is slow in this material. A possible solution to increase the ion diffusion is by reducing the Ca–Ca distance in this material via synthesizing nanoparticles.⁴⁶ Table II lists the activation energies calculated for Ca-ion migration in different Ca-based oxide materials.^12,13,47 Torres et al.¹³ recently performed DFT calculations on a variety of Ca based minerals and reported the activation energies of Ca-ion migration. Interestingly, Ca₃Cr₂Si₃O₁₂ and Ca₃Mn₂Si₃O₁₂ materials, isostructural with Ca₃Fe₂Si₃O₁₂, were also considered, and their activation energies of migration were reported to be 2.07 eV and 2.09 eV, respectively.¹³ Activation energies calculated for Ca-ion migration are lower by ∼0.50 eV than those calculated in a similar Ca₃Fe₂Si₃O₁₂ structure. This is due to different methodologies considered and different transition elements present in the crystal structures. Activation energies are greater than 2 eV in all cases, inferring poor ionic conductivity. A possible reason can be due to the strong electrostatic interaction between the migrating Ca²⁺ cations and the other ions in the lattice.

This section examines the introduction of cations via divalent, trivalent, and tetravalent doping into Ca₃Fe₂Si₃O₁₂. Doping can result in beneficial or detrimental effects on the material. In particular, suitable dopants can stabilize a material in order to prevent unfavorable phase transformations. The aim here is to provide guidance to the dopants that require the least energy to integrate into the lattice, the effects of which can be studied in subsequent experimental testing.

Divalent cations (M = Ni, Mg, Co, Mn, Sr, and Ba) were first considered for doping on the Ca site. The doping mechanism can be seen in the following equation:

The calculated solution enthalpies are reported in Fig. 4. The lowest solution enthalpy (0.21 eV) is predicted for Mn²⁺, inferring the highest solubility. The lowest solution enthalpy suggests that Mn²⁺ would be an ideal candidate to stabilize Ca₃Fe₂Si₃O₁₂. The ionic radius of Ca²⁺ in an eight coordination environment is 1.00 Å. Favorability of the Mn²⁺ dopant can be due to its ionic radius (0.96 Å), which is close to the ionic radius of the Ca²⁺ ion. Both Mg²⁺ and Co²⁺ exhibit an identical solution enthalpy of 0.59 eV. This is due to a very small difference in their ionic radii. Solution enthalpy increases with the difference in the ionic radius between dopants and the Ca²⁺ ion. The solution enthalpies of Sr²⁺ and Ni²⁺ were calculated to be 0.64 eV and 0.77 eV, respectively. The highest solution enthalpy of 2.20 eV is calculated for Ba²⁺, suggesting that this dopant is very unfavorable for doping on the Ca site.

To increase the capacity of Ca₃Fe₂Si₃O₁₂, divalent dopants were considered on the Fe site. The following reaction equation [Eq. (16)] shows that this doping strategy can introduce additional Ca²⁺ ions in the form of Ca interstitials in the lattice. A previous experimental study considered Co³⁺ doping on the Ru site in Li₂RuO₃ to increase the Li⁺ concentration, and doping resulted in an enhancement in the electrochemical reversibility of Li⁺ ions and the concentration of Li⁺ ions in the lattice,⁴⁸ 

The favorable dopant for this process is Mn²⁺, and its solution energy is calculated to be 3.36 eV. Solution energies calculated for Co²⁺, Ni²⁺, Mg²⁺, and Zn²⁺ are very close to the value calculated for Mn²⁺; therefore, they are also candidate dopants for experimental testing (refer to Fig. 5). The ionic radius of Fe³⁺ is 0.65 Å. Favorability of the five dopants (Mn, Co, Ni, Mg, and Zn) is due to their ionic radii (0.65 Å–0.74 Å) matching closely with the ionic radius of Fe³⁺. The ionic radius of Sr²⁺ (1.18 Å) deviates significantly from that of Fe³⁺, reflecting in solution enthalpy. The highest positive solution enthalpy of 6.17 eV is calculated for Ba²⁺, suggesting that this dopant is highly unfavorable on the Fe site.

Here, we consider a range of trivalent cations (M = Al, Ga, In, Sc, Y, Gd, and La) on the Fe and Si sites. Doping on the Fe site produced no charge compensating defects as the charge on Fe is +3 [refer to Eq. (17)]. Solution enthalpies are plotted as a function of ionic radii of dopants in Fig. 4(b),

The most favorable dopant is Sc³⁺, with a solution enthalpy of 0.09 eV. The solution enthalpy of In³⁺ is higher only by 0.02 eV than that calculated for Sc³⁺, meaning that In³⁺ is also a promising dopant. Favorability of these two dopants can be partly due to their ionic radii being closer to the ionic radius of Fe³⁺ (0.65 Å). The solution enthalpy of Ga³⁺ is 0.29 eV although its ionic radius (0.62 Å) is very close to the ionic radius of Fe³⁺. Other dopants exhibit high solution enthalpies due to their ionic radii deviating from the ionic radius of Fe³⁺ [refer to Fig. 4(b)].

The formation of additional Ca²⁺ ions in the lattice was considered by doping trivalent cations on the Si site. As we discussed earlier, this process can also introduce Ca interstitials, as shown in the following reaction,

The solution enthalpies are reported in Fig. 5(b). In all cases, solution enthalpies are highly positive (>5 eV), suggesting that this process required high energy possibly in the form of heat. The most favorable dopant is Al³⁺, with a solution enthalpy of 5.49 eV. The lowest solution enthalpy for Al³⁺ can be due to its ionic radius (0.39 Å) being closer to the ionic radius of Si⁴⁺ (0.26 Å). High endoergic solution enthalpy for this mechanism can be due to the high defect formation energy of quadruply charged Si.

Finally, tetravalent dopants (M = Ge, Ti, Sn, Zr, and Ce) were considered on the Si site. Solution enthalpies are reported in Fig. 4(c). The following defect equation explains the mechanism involved in this process,

The most favorable solution energy (0.43 eV) is calculated for Ge⁴⁺. This is due to the ionic radius of Ge⁴⁺ (0.39 Å) being closer to that of Si⁴⁺ (0.26 Å). The solution enthalpy for Ti⁴⁺ is 3.86 eV, meaning that it is an unfavorable dopant. The second favorable dopant is Sn⁴⁺, with a solution enthalpy of 2.10 eV. Solution enthalpy gradually increases with an increase in the ionic radius. The highest enthalpy of solution for Ce⁴⁺ (4.12 eV) suggests that doping would require external heat energy.

Here, we calculate the activation energies for Ca-ion diffusion in the presence of Mn²⁺ on the Fe site and Al³⁺ on the Si site. Figure 6 shows the energy profile diagrams for the local Ca hops. There is a reduction (by 0.45 eV) in the activation energy for the doping of Mn²⁺ on the Fe site [refer to Fig. 6(a)]. In the case of doping of Al³⁺ on the Si site, activation energy is reduced only by 0.09 eV [refer to Fig. 6(b)]. In both cases, Ca–Ca distances have been slightly elongated compared to the distances present in the un-doped crystal structure. This perturbation in distances is due to the charge and ionic radius mismatch between Mn²⁺ and Fe³⁺ and Al³⁺ and Si⁴⁺. The current simulation shows that doping of Mn²⁺ and Al³⁺ would simultaneously increase the concentration of Ca²⁺ ions in the lattice and reduce the activation energy of the Ca-ion migration.

DFT simulations were used to examine the electronic structures of doped configurations. Here, we discuss the results of the most favorable dopants (Mn²⁺, Ga³⁺, Al³⁺, and Ge⁴⁺), as discussed in the previous sections.

The DOS plot calculated for bulk Ca₃Fe₂Si₃O₁₂ is shown in Fig. 6(a). Metallic behavior is observed for the ferromagnetic configuration of bulk Ca₃Fe₂Si₃O₁₂. Doping of Mn alters the electronic structure slightly, as shown in the DOS plot [refer to Fig. 7(b)]. Peaks arising from the 3d states of Mn appear near the Fermi level [refer to Fig. 7(c)]. Calculated Mn–O bond distances (2.20 Å–2.50 Å) [refer to Fig. 7(d)] are slightly shorter than the Ca–O (2.37 Å–2.54 Å) bond distances due to the ionic radius of Mn²⁺ being smaller than that of Ca²⁺. The constant charge density plot of the Mn-doped configuration is shown in Fig. 7(e).

In the case of Sc, the Fermi level does not change significantly, and the doped configuration maintains its metallic behavior [refer to Fig. 7(f)]. The Fermi level does not have any contribution from Sc. States associated with Sc appear at ∼8 eV in the conduction band [refer to Fig. 7(g)]. Sc–O bond distances are slightly longer (by ∼0.04 Å) than the Fe–O bond distances. This is due to the ionic radius of Sc³⁺ being larger than that of Fe³⁺, as discussed in Sec. III D 2.

Doping of Ge on the Si site has a little effect on the DOS structure [refer to Fig. 7(j)]. The doped configuration exhibits metallic behavior. Small peaks that belong to the 3s and 3p states of Ge appear just below the Fermi level [refer to Fig. 7(k)]. The Si–O bond distance in pristine Ca₃Fe₂Si₃O₁₂ is 1.66 Å. The Ge–O bond distances are longer than (by 0.15 Å) the Si–O bond distances. This is because of the ionic radius of Ge⁴⁺ (0.39 Å) being larger than that of Si⁴⁺ (0.26 Å).

Doping of Al on the Si site shifts the Fermi energy level only by 0.02 eV [refer to Fig. 8(b)]. The doped configuration is metallic. Additional peaks arising from 3s and 3p states of Al appear near the Fermi level [refer to Fig. 8(c)]. The calculated Al–O bond distances (1.77 Å) are longer than the Si–O bond distances (1.66 Å). This is because of the ionic radius and charge density of Al³⁺ being larger and smaller, respectively, than that of Si⁴⁺.

There is a slight shift in the Fermi energy level by 0.14 eV upon doping of Mn on the Fe site. Peaks associated with Mn appear near the Fermi level and in the valence band. Metallic behavior is still maintained. A very small perturbation is observed in the Mn–O bond distances compared to Fe–O bond distances. This is due to the small difference in the ionic radius and charge mismatch. The constant charge density plot of Mn doped on the Fe site in Ca₃Fe₂Si₃O₁₂ is shown in Fig. 8(h).

Computational simulation was applied to analyze the behavior of the defects, diffusion pathways, and dopant properties in Ca₃Fe₂Si₃O₁₂. The Ca/Fe anti-site defect was predicted to be the most energetically favorable intrinsic defect process. The activation energy for three dimensional Ca vacancy migration in Ca₃Fe₂Si₃O₁₂ was calculated to be 2.63 eV, inferring low calcium mobility. Low solution energies were calculated for Mn²⁺, Sc³⁺, and Ge⁴⁺ on the Ca, Fe, and Si sites, respectively. Here, we show that Al³⁺ and Mn²⁺ are promising dopants on the Si and Fe sites, respectively, to introduce additional Ca in the form of Ca interstitials in Ca₃Fe₂Si₃O₁₂ for improving its capacity and Ca-ion diffusion. The difference in the electronic properties between un-doped and doped configurations was explained using electronic structures calculated using DFT.

See the supplementary material for the Buckingham potentials used in the classical simulation of the current study.

The data that support the findings of this study are available from the corresponding author upon reasonable request.

Computational facilities and support were provided by the High Performance Computing Centre at Imperial College London.

This research was financially supported by the European Union’s H2020 Programme under Grant Agreement No. 824072–HARVESTORE.

The authors declare that there is no competing financial interest.