Defects and Calcium Di ﬀ usion in Wollastonite

: Wollastonite (CaSiO 3 ) is an important mineral that is widely used in ceramics and polymer industries. Defect energetics, di ﬀ usion of Ca ions and a solution of dopants are studied using atomistic-scale simulation based on the classical pair potentials. The energetically favourable defect process is calculated to be the Ca-Si anti-site defect cluster in which both Ca and Si swap their atomic positions simultaneously. It is calculated that the Ca ion migrates in the ab plane with an activation energy of 1.59 eV, inferring its slow di ﬀ usion. Favourable isovalent dopants on the Ca and Si sites are Sr 2 + and Ge 4 + , respectively. Subvalent doping by Al on the Si site is a favourable process to incorporate additional Ca in the form of interstitials in CaSiO 3 . This engineering strategy would increase the capacity of this material.


Introduction
Wollastonite is an important naturally occurring or synthetic mineral of great interest in the development of ceramics, plastics and paints [1][2][3][4]. This mineral is mainly found in the USA, India, Mexico, China, Turkey and Finland [5]. In general, a relatively small amount of impurities, such as Fe, Mg, Al and Sr, contaminate pure CaSiO 3 [6]. Owing to its remarkable physical, mechanical, electrical and thermal properties, this material has been used as a ceramic, an insulator, a filler and a polymer [2,[7][8][9][10][11][12].
Though wollastonite is mined from ores for large-scale commercial applications, there are many experimental reports explaining the procedures of synthesising and characterising wollastonite in laboratories [1,5,10,13,14]. Abd Rashid et al. [1] used the solid-state reaction method at low temperature to produce wollastonite from limestone and silica sand and concluded that the resultant compound is expected to exhibit good physical properties due to the experimental density that is close to its theoretical value. The precipitation technique is the most common method that has been reported by several researchers to prepare wollastonite [15][16][17][18]. The sol-gel method is another technique that has been used to synthesise wollastonite by De La Casa-Lillo et al. [19], and the correlation between thermal treatment and bio activity has been discussed. In the literature, a few theoretical studies are available on wollastonite [20][21][22][23]. Edrees et al. [20] used density functional theory (DFT) simulations to analyse the structural, mechanical, optical and electronic properties of wollastonite. Atomistic simulations performed on the surface structures of wollastonite show that the stabilisation of the surface is due to the adsorption of water in a dissociated form [21]. Longo et al. [22] used quantum mechanical simulations to show that the adsorption of CO 2 on the CaSiO 3 surface facilitates the formation of CO 3 2− ions due to the reaction between CO 2 and surface oxygen. DFT simulations performed by Profeta et al. [23] show that accurate 17 O NMR spectra should be calculated using hybrid functions. However, there are no experimental or theoretical studies on intrinsic defects, diffusion or dopants of wollastonite in the literature. The presence of defects in a material is important, as they influence its physical, mechanical and optical properties. Atomistic simulation based on the pair-wise potentials is a powerful tool to examine the energetics of intrinsic defects, migration pathways and solutions of dopants in ionic materials. A wide range of oxide materials have been examined using this method to provide insight into defects, diffusion and dopants [24][25][26][27][28]. The present study systematically examined the intrinsic defect processes, Ca-ion diffusion and solution of divalent (Co, Mn, Ni, Mg, Zn, Sr and Ba), trivalent (Al, Ga In, Fe, Sc, Y, Gd and La) and tetravalent (Ga, Ti, Sn, Zr and Ce) dopants in CaSiO 3 .

Computational Methods
A classical pair potential-based atomistic simulation, as implemented in the General Utility Lattice Program (GULP) code [29], was used to calculate defect, diffusion and dopant properties in CaSiO 3 . Interactions between ions were modelled using long-range (Coulomb) and short-range (electron-electron repulsion and dispersive attraction) forces. Short-range forces were modelled using Buckingham potentials (refer to Table 1). Full geometry optimisation (cell parameters and ionic positions) was performed using the Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm [30]. The Mott-Littleton approach was implemented to model the defects [31]. The diffusion of Ca ions was calculated considering two adjacent Ca vacancy sites as initial and final configurations. Seven atomic positions were selected between two adjacent Ca vacancy sites and they were allowed to relax. The activation energy is defined as the difference between the vacancy formation energy and the maximum energy along the diffusion path. In the current methodology, ions were fully charged and defect concentration was low. Nevertheless, relative energies and trends should be consistent [32]. Table 1. Buckingham potential parameters [33,34] used in the classical simulations of CaSiO 3 . Two-body , where A, ρ and C are parameters which were selected carefully to reproduce the experimental data. The values of Y and K represent the shell charges and spring constants. A very large spring constant means there is no shell charge and atom is treated as core.

Crystal Structure of CaSiO 3 and Validation of Potentials
The crystal structure of CaSiO 3 is monoclinic [space group P2 1 /a (no: 14)] with the experimental lattice parameters a = 15.41 Å, b = 7.32 Å, c = 7.06 Å, α = γ = 90 • and β = 95.3 • [35]. Ca ions form a distorted octahedral structure and Si atoms form tetrahedral units. SiO 4 tetrahedral units are inter-connected by corner sharing (refer to Figure 1). Geometry optimisation was performed under constant pressure to get equilibrium lattice constants and compare them with experimental values to validate the potentials used in this study. Calculated lattice parameters are reported in Table 2. There was a good agreement between calculated and experimental values.

Defect Energetics
In this section, we discuss the energetics of key defect processes in CaSiO3. Point defect (vacancies and interstitials) energies were first calculated and then they were combined to calculate Schottky and Frenkel energies. Anti-site defects in which Ca and Si swap their atomic positions were also calculated in the form of isolated and cluster. The following equations in Kröger-Vink notation [36] describe the Schottky, Frenkel and anti-site defect process.
Ca/Si antisite (cluster): Ca + Si → {Ca : Si •• } X (8) Figure 2 reports the defect reaction energies. The Frenkel energies are important, as they influence the diffusion properties of the material. Furthermore, a low Frenkel energy is associated with the high concentration of point defects, leading to the loss of the crystal structure. The oxygen Frenkel energy is the lowest energy (2.33 eV) process among the other Frenkel energies. The Si Frenkel energy (8.95 eV) is significantly higher than the Ca Frenkel energy (3.46 eV). This is due to the introduction of highly charged defects ( and Si •••• ) in the lattice. The formation of Schottky

Defect Energetics
In this section, we discuss the energetics of key defect processes in CaSiO 3 . Point defect (vacancies and interstitials) energies were first calculated and then they were combined to calculate Schottky and Frenkel energies. Anti-site defects in which Ca and Si swap their atomic positions were also calculated in the form of isolated and cluster. The following equations in Kröger-Vink notation [36] describe the Schottky, Frenkel and anti-site defect process.
Ca/Si antisite (cluster) : Chemistry 2020, 2 940 Figure 2 reports the defect reaction energies. The Frenkel energies are important, as they influence the diffusion properties of the material. Furthermore, a low Frenkel energy is associated with the high concentration of point defects, leading to the loss of the crystal structure. The oxygen Frenkel energy is the lowest energy (2.33 eV) process among the other Frenkel energies. The Si Frenkel energy (8.95 eV) is significantly higher than the Ca Frenkel energy (3.46 eV). This is due to the introduction of highly charged defects (V Si and Si •••• i ) in the lattice. The formation of Schottky defects is endoergic by 3.58 eV, meaning that this defect is unfavourable. The formation of CaO and SiO 2 via CaO Schottky and SiO 2 Schottky processes, respectively, was also considered. Their defect energies are also highly endothermic. Finally, the Ca/Si anti-site defect was considered. The cluster form of this defect exhibits lower energy than its isolated form. This is because of the unstable nature of isolated defects (Ca Si and Si •• Ca ) aggregating to form clusters with a binding energy of −2.53 eV. The anti-site defect has been observed in many oxide materials during the cycling of as-prepared materials and synthesis at high temperatures and pressures [37][38][39]. This defect has also been observed in previous theoretical studies [40][41][42][43].
Chemistry 2020, 2, x 4 defects is endoergic by 3.58 eV, meaning that this defect is unfavourable. The formation of CaO and SiO2 via CaO Schottky and SiO2 Schottky processes, respectively, was also considered. Their defect energies are also highly endothermic. Finally, the Ca/Si anti-site defect was considered. The cluster form of this defect exhibits lower energy than its isolated form. This is because of the unstable nature of isolated defects (Ca and Si •• ) aggregating to form clusters with a binding energy of -2.53 eV. The anti-site defect has been observed in many oxide materials during the cycling of as-prepared materials and synthesis at high temperatures and pressures [37][38][39]. This defect has also been observed in previous theoretical studies [40][41][42][43].

Calcium Ion Diffusion
The performance of a material can be influenced by its diffusion properties, which are associated with the activation energy and pre-exponential factor of the Arrhenius equation and energetics of point defects. The present methodology is an appropriate tool to calculate the activation energies of long-range diffusion pathways of Ca ions. The results presented here would be beneficial to the experimentalist as the determination of ion pathways and activation energies can be challenging experimentally. In previous modelling studies with this classical approach, many oxides have been examined and diffusion properties reported [27,[44][45][46]. For example, the diffusion pathway calculated in LiFePO4 by Fisher et al. [47] is in excellent agreement with the neutron diffraction study reported by Nishimura et al. [48].
In general, Ca ion diffusion in Ca-bearing oxide materials is expected to be sluggish due to its ions with double positive charge. In previous theoretical simulations based on the DFT and classical pair potentials, it was shown that activation energies of Ca ion diffusion in Ca-based minerals are high, confirming the slow diffusion [49][50][51][52].
We identified four different local Ca hops (A, B, C and D) and their vacancy-assisted migration pathways were calculated (refer to Figure 3). Table 3 reports the local Ca hops with their distances and corresponding activation energies.

Calcium Ion Diffusion
The performance of a material can be influenced by its diffusion properties, which are associated with the activation energy and pre-exponential factor of the Arrhenius equation and energetics of point defects. The present methodology is an appropriate tool to calculate the activation energies of long-range diffusion pathways of Ca ions. The results presented here would be beneficial to the experimentalist as the determination of ion pathways and activation energies can be challenging experimentally. In previous modelling studies with this classical approach, many oxides have been examined and diffusion properties reported [27,[44][45][46]. For example, the diffusion pathway calculated in LiFePO 4 by Fisher et al. [47] is in excellent agreement with the neutron diffraction study reported by Nishimura et al. [48].
In general, Ca ion diffusion in Ca-bearing oxide materials is expected to be sluggish due to its ions with double positive charge. In previous theoretical simulations based on the DFT and classical pair potentials, it was shown that activation energies of Ca ion diffusion in Ca-based minerals are high, confirming the slow diffusion [49][50][51][52].
We identified four different local Ca hops (A, B, C and D) and their vacancy-assisted migration pathways were calculated (refer to Figure 3). Table 3 reports the local Ca hops with their distances and corresponding activation energies. Chemistry 2020, 2, x 5 The local hop A exhibits the lowest activation energy. The activation energy calculated for the local hop B is slightly higher by 0.03 eV than that calculated for the local hop A. The highest activation energy is noted for the local hop D. In order to find the long-range diffusion pathways, local hops were connected. Three possible long-range diffusion pathways were identified (refer to Table 4). The lowest energy long-range diffusion pathway is A↔B↔C↔A↔B with an activation energy of 1.59 eV, inferring the slow diffusion. The movement of Ca ions is observed in the ab plane. Figure 4 shows the energy profile diagrams for local hops.   The local hop A exhibits the lowest activation energy. The activation energy calculated for the local hop B is slightly higher by 0.03 eV than that calculated for the local hop A. The highest activation energy is noted for the local hop D. In order to find the long-range diffusion pathways, local hops were connected. Three possible long-range diffusion pathways were identified (refer to Table 4). The lowest energy long-range diffusion pathway is A↔B↔C↔A↔B with an activation energy of 1.59 eV, inferring the slow diffusion. The movement of Ca ions is observed in the ab plane. Figure 4 shows the energy profile diagrams for local hops.

Solution of Dopants
Dopants play a significant role in tailoring the properties of a material. In particular, dopants with different sizes or charges compared to the host atoms will affect the properties of the host compound. Here, we considered a variety of isovalent and aliovalent dopants for screening and predicting promising dopants that can be verified and further explored experimentally. Solution energies were calculated using appropriate charge compensation defects and lattice energies. The Supplementary Materials provide the Buckingham potentials used for dopants (refer to Table S1).

Divalent Dopants
First, the Ca site was considered to dope divalent dopants (M = Co, Mn, Ni, Mg, Zn, Sr and Ba). The solution energy was calculated using the following equation.
Exothermic solution energy (−0.03 eV) was calculated for Sr 2+ , suggesting that it is the most favourable dopant on the Ca site (refer to Figure 5). A possible composition that can be prepared by experiment is Ca 1-x Sr x SiO 3 (0.0 < x < 1.0). The second most stable dopant is Ba 2+ , with a solution energy of 0.16 eV. Other dopants exhibit high endoergic solution energies, meaning that they are unlikely to be doped at room temperature. The high preference of Sr 2+ could be due to its ionic radius (1.18 Å), which is close to the ionic radius of Ca 2+ (1.00 Å).

Solution of Dopants
Dopants play a significant role in tailoring the properties of a material. In particular, dopants with different sizes or charges compared to the host atoms will affect the properties of the host compound. Here, we considered a variety of isovalent and aliovalent dopants for screening and predicting promising dopants that can be verified and further explored experimentally. Solution energies were calculated using appropriate charge compensation defects and lattice energies. The Supplementary Materials provide the Buckingham potentials used for dopants (refer to Table S1).

Divalent Dopants
First, the Ca site was considered to dope divalent dopants (M = Co, Mn, Ni, Mg, Zn, Sr and Ba). The solution energy was calculated using the following equation.

MO + Ca → M + CaO
Exothermic solution energy (-0.03 eV) was calculated for Sr 2+ , suggesting that it is the most favourable dopant on the Ca site (refer to Figure 5). A possible composition that can be prepared by experiment is Ca1-xSrxSiO3 (0.0 < x < 1.0). The second most stable dopant is Ba 2+ , with a solution energy of 0.16 eV. Other dopants exhibit high endoergic solution energies, meaning that they are unlikely to be doped at room temperature. The high preference of Sr 2+ could be due to its ionic radius (1.18 Å), which is close to the ionic radius of Ca 2+ (1.00 Å).

Trivalent Dopants
The doping of trivalent cations (Al 3+ , Ga 3+ , In 3+ , Fe 3+ , Sc 3+, Y 3+ , Gd 3+ and La 3+ ) was considered on the Si site. In this doping process, Ca interstitials were introduced to compensate the negatively charged lattice. This process can increase the capacity of CaSiO3 and enhance the diffusion of Ca ions.

Trivalent Dopants
The doping of trivalent cations (Al 3+ , Ga 3+ , In 3+ , Fe 3+ , Sc 3+, Y 3+ , Gd 3+ and La 3+ ) was considered on the Si site. In this doping process, Ca interstitials were introduced to compensate the negatively charged lattice. This process can increase the capacity of CaSiO 3 and enhance the diffusion of Ca ions. The following equation was used to calculate the solution energy for this process.
Calculated solution enthalpies are reported in Figure 6. The most favourable dopant for this process is Al 3+ , with a solution energy of 3.47 eV, though the solution is endoergic. The possible composition that can be prepared by experiment would be Ca 1+x Al x Si 1−x O 8 (0.0 < x< 1.0). The second most favourable dopant is Fe 3+ . High solution energies (>4 eV) are calculated for the other dopants, meaning they are unlikely to take place at room temperature.

MO + Si → M + SiO
Solution energies are reported in Figure 7. The promising dopant for this process is Ge 4+ . The preference of Ge 4+ is due to the ionic radius of Si 4+ (0.26 Å), which is close to the ionic radius of Ge 4+ (0.39 Å). Endothermic solution energy shows that this process requires energy in the form of heat. This is due to the stronger Si-O bonds as compared to the Ge-O bonds. The second most favourable dopant is Ti 4+ . Its solution energy is higher only by 0.07 eV than that calculated for Ge 4+ . High solution energies are calculated for the other dopants. In particular, solution energies for ZrO2 and CeO2 are 2.37 eV and 3.30 eV, respectively, suggesting that they are unlikely to occur at room temperature.
Solution energies are reported in Figure 7. The promising dopant for this process is Ge 4+ . The preference of Ge 4+ is due to the ionic radius of Si 4+ (0.26 Å), which is close to the ionic radius of Ge 4+ (0.39 Å). Endothermic solution energy shows that this process requires energy in the form of heat. This is due to the stronger Si-O bonds as compared to the Ge-O bonds. The second most favourable dopant is Ti 4+ . Its solution energy is higher only by 0.07 eV than that calculated for Ge 4+ . High solution energies are calculated for the other dopants. In particular, solution energies for ZrO 2 and CeO 2 are 2.37 eV and 3.30 eV, respectively, suggesting that they are unlikely to occur at room temperature.
Solution energies are reported in Figure 7. The promising dopant for this process is Ge 4+ . The preference of Ge 4+ is due to the ionic radius of Si 4+ (0.26 Å), which is close to the ionic radius of Ge 4+ (0.39 Å). Endothermic solution energy shows that this process requires energy in the form of heat. This is due to the stronger Si-O bonds as compared to the Ge-O bonds. The second most favourable dopant is Ti 4+ . Its solution energy is higher only by 0.07 eV than that calculated for Ge 4+ . High solution energies are calculated for the other dopants. In particular, solution energies for ZrO2 and CeO2 are 2.37 eV and 3.30 eV, respectively, suggesting that they are unlikely to occur at room temperature.

Conclusions
In this study, an atomistic simulation based on the classical potentials was applied to examine the intrinsic defect processes, vacancy-assisted Ca ion diffusion and solution of dopants in CaSiO3. The Ca-Si anti-site defect cluster was calculated to be the lowest energy process, meaning that a small amount of Ca on the Si site and Si on the Ca site will be present, particularly at higher temperatures. The long-range vacancy-assisted Ca ion migration pathway is in the ab plane with an activation energy of 1.59 eV, inferring slow diffusion. The promising dopant on the Si site to increase the Ca

Conclusions
In this study, an atomistic simulation based on the classical potentials was applied to examine the intrinsic defect processes, vacancy-assisted Ca ion diffusion and solution of dopants in CaSiO 3 . The Ca-Si anti-site defect cluster was calculated to be the lowest energy process, meaning that a small amount of Ca on the Si site and Si on the Ca site will be present, particularly at higher temperatures. The long-range vacancy-assisted Ca ion migration pathway is in the ab plane with an activation energy of 1.59 eV, inferring slow diffusion. The promising dopant on the Si site to increase the Ca content in CaSiO 3 is Al 3+ . Promising isovalent dopants on the Ca and Si sites are Sr 2+ and Ge 4+ , respectively, and these dopants can prevent phase transformation.