Defects, Di ﬀ usion and Dopants in Sillimanite

: Aluminum silicate based mineral “Sillimanite” (Al 2 SiO 5 ) is important in the industrial preparation of aluminum-silicon alloys and cement. In the present study classical pair potential simulations are used to examine the intrinsic defect processes, di ﬀ usion pathways of Al 3 + and O 2 − ions together with their activation energies and promising dopants on the Al and Si sites in Al 2 SiO 5 . The cation anti-site (Al-Si) defect cluster is calculated to be the most favorable defect, highlighting the cation disorder in this material, in agreement with the experiment. The cation disorder is important as this defect can change the mechanical and chemical properties of Al 2 SiO 5 . The Al 3 + ions and O 2 − ions migrate in the c direction with corresponding activation energies of 2.26 eV and 2.75 eV inferring slow ion di ﬀ usion. The prominent isovalent dopants on the Al and Si sites are found to be the Ga and Ge, respectively, suggesting that they can be used to prevent phase transformation and tune the properties of sillimanite.


Introduction
Sillimanite (Al 2 SiO 5 ) is a naturally occurring mineral found in many parts of the world [1][2][3][4][5]. This mineral is mainly found in metamorphosed rocks and its formation is dependent on temperature and pressure [6]. It is an economically important mineral in the industrial preparation of bricks, cement, ceramics, jewelry (e.g., sillimanite gold ring) and fine porcelain (e.g., table top) [7][8][9]. The industrial use of sillimanite is related to its unique chemical composition, thermal stability and the formation of mullite-rich aggregates [10].
In order to produce high quality sillimanite samples from low grade ores to meet industrial needs, a variety of techniques have been applied [4,10]. Jin et al. [5] used a flotation technique to effectively separate sillimanite from ores and concluded that a higher flotation recovery is observed for sillimanite than its polymorph "kyanite" in the presence of a sodium oleate collector. The composition of sillimanite from high-grade metamorphic rock was analyzed by Grew et al. [11] and it was shown that a small amount of Fe 3+ ions are present on the Al site forming stoichiometric (Al,Fe 3+ ) 2 SiO 5 composition. Annealing experiments undertaken by Holland et al. [12] at high pressures, resulted an Al/Si disorder in sillimanite.
Though there are many experimental studies on natural sillimanite, only a few experimental studies are available on the synthetic sillimanite. Xu et al. [13] synthesized sillimanite in the form of whiskers at low temperature and characteried its structure using X-ray diffraction (XRD) together with electron microscopy techniques. One of the difficulties in the preparation of sillimanite is the phase transformation in which sillimanite crystal structure becomes distorted and in this distorted structure (mullite), a disordered distribution of Si and Al is present. Igami et al. [14] used synchrotron X-ray diffraction experiments to determine the temperature at which sillimanite transformation occurs. It was concluded that the mullitization (sillimanite to mullite) temperature is at~1200 • C [14]. As synthetic sillimanite is of interest in the ceramic industry, its defect properties are also important in order to optimize its properties. Defect studies on this material have not been explored by experiments yet. Theoretical studies can provide useful information on the defects in this material. To the best of our knowledge, there are no theoretical studies available in the literature on the defects, including cation disorder (anti-site defect) and dopants in sillimanite.
Atomistic simulation techniques based on the classical pair potentials can provide detailed information about defect energetics, cation disorder and diffusion pathways, together with activation energies and promising dopants. Such information is useful in the interpretation of experimental data. In particular, the cation disorder and substitutional doping are important as these defect processes can influence the mechanical and chemical properties of Al 2 SiO 5 . The current methodology has been successfully applied to a variety of oxide materials to make precise predictions on the defect structures, local structural changes, migration energies and dopant properties [15][16][17][18][19][20]. In this work, atomistic simulation techniques are used to examine various defect processes, diffusion pathways and isovalent dopant behavior in sillimanite.

Computational Methods
All calculations were performed using a classical simulation code General Utility Lattice Program (GULP, version 3.4.1) [21]. Ionic interactions were described using long range (Coulombic) and short range (Pauli repulsion and van der Waals attraction) forces. The short range interactions were modelled using Buckingham potentials. Energy minimization calculations, to relax the simulation boxes and ionic positions, were performed using the Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm [22], as implemented in the GULP code. A gradient norm of 0.001 eV/Å was used to converge bulk and defect calculations. The Mott-Littleton method [23] was employed to calculate defect energies. In this approach, the crystal lattice is divided into two spherical regions (region I and region II). The inner spherical region (region I) consisted of~700 ions immediately surrounding the defect and this region was relaxed explicitly. The ions in the outer sphere (region II) were relaxed using quasi-continuum methods. Two adjacent Al vacancy sites were first created, and Al-ion interstitial positions were then systematically placed at regular intervals along the diagonal connecting them. Seven interstitial positions were considered in all cases and the interstitial ion was fixed while all other ions were free to relax. However, fixing the interstitial ion position does not guarantee the minimum energy path, and it will give only a direct diffusion path. Therefore, interstitial positions were allowed to move in the x, y, z, xy, yz, and xz directions separately. Finally, the lowest activation energy pathway (curved pathway) was reported. The activation energy of the diffusing ion was calculated by taking the energy difference between the initial configuration and the saddle point configuration. This methodology has been discussed in previous theoretical studies [24,25].

Crystal Structure of Al 2 SiO 5
Sillimanite exhibits orthorhombic structure (space group Pnma (no 62)) with unit cell parameters a = 7.4856, b = 7.6738, c = 5.7698 Å and α = β = γ = 90.0 • [26]. While tetrahedral coordination is observed for Si 4+ ions, both tetrahedral and octahedral coordination are noted for Al 3+ ions. Chains of edge-sharing AlO 6 octahedra exist in this structure along the (001) direction (refer to Figure 1). In the same direction, alternating tetrahedral AlO 4 and SiO 4 units are present, sharing their corners. Buckingham potentials (refer to Table 1) [27,28] were validated by comparing experimental lattice parameters with those calculated from geometry optimization of bulk Al 2 SiO 5 . There is a good agreement between calculated and experimental values, as reported in Table 2.   [27,28] used in the classical simulations of Al2SiO5. Twobody Φij (rij) = Aij exp (−rij/ρij) − Cij/rij 6 , 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. The main diagonal components of the elastic tensors were also calculated and compared with available theoretical and experimental values to further validate the potential parameters. Table 3 lists the calculated and experimental values. In general, there is a reasonable agreement between calculated values in this study and the values from the other theoretical studies and the experiment. In order to validate the potentials against energetics, we performed energy minimization calculations on the two other polymorphs of Sillimanite (Andalusite and Kyanite). Relative energies  Table 1. Buckingham potential parameters [27,28] used in the classical simulations of Al 2 SiO 5 . Two-body The main diagonal components of the elastic tensors were also calculated and compared with available theoretical and experimental values to further validate the potential parameters. Table 3 lists the calculated and experimental values. In general, there is a reasonable agreement between calculated values in this study and the values from the other theoretical studies and the experiment. In order to validate the potentials against energetics, we performed energy minimization calculations on the two other polymorphs of Sillimanite (Andalusite and Kyanite). Relative energies were compared with the experimental study reported by Waldbaum [32]. There is a good agreement in the trend between the calculation and the experiment (refer to Table 4).

Intrinsic Defect Energetics
In this section, a systematic survey of intrinsic defects in Al 2 SiO 5 was performed. A series of point defects (vacancies and interstitials) were first calculated and then they were combined to calculate Schottky and Frenkel energies. The following equations (Equations (1)-(8)) as written using Kröger-Vink notation [33] describe the Schottky, Frenkel and anti-site defect processes.
Al/Si antisite (cluster) : Calculations show that the lowest energy defect process is the Al-Si anti-site defect cluster (0.43 eV/defect) (refer to Figure 2). This indicates that a small concentration of cation mixing (Al-Si) is present in Al 2 SiO 5 . This is in agreement with the experimental study by Holland et al. [12]. The defect energy for the isolated form of the anti-site defect is 0.82 eV. Therefore, the energy difference between the isolated and cluster forms of the anti-site defect is −0.39 eV. This energy is the binding energy and the exoergic nature of the binding energy means that the isolated anti-site defects are not stable and they aggregate to form cluster without energy penalty. The anti-site defect has been identified in a variety of oxide materials experimentally and theoretically [34][35][36][37][38]. For example, Li/Fe anti-site disorder was determined experimentally during the cycling of Li 2 FeSiO 4 cathode material for Li-ion batteries [39]. The next favorable defect is the O-Frenkel (1.98 eV/defect). This defect would not be present at a significant concentration at normal temperatures. The Al 2 O 3 Schottky defect energy is 2.97 eV/defect, suggesting that loss of Al 2 O 3 is possible only at high temperatures. The other Frenkel and Schottky defects exhibit high formation energies and they would not be present under any conditions.

Self-Diffusion of Aluminium and Oxygen
The ionic conductivity of a material is mainly dependent on the parameters associated with the formation of point defects, pre-exponential factor and activation energy of the migrating ion. The current simulation technique enables the examination of various ion diffusion pathways with their activation energies in Al2SiO5. Experimental determination of ion diffusion pathways is often difficult. The current methodology has successfully reproduced experimentally determined ion diffusion pathways in a variety of ionic oxide materials [16]. For example, an excellent agreement between the calculated and experimental Li-ion diffusion pathways was observed in LiFePO4 [40,41].
Vacancy assisted diffusion pathways were examined for Al 3+ and O 2− ions as their Frenkel energies are lower than the other Frenkels. Two different local Al hops (A and B) (refer to Figure 3) with the jump distances of 2.97 Å and 3.14 Å were identified (refer to Figure 3). Individual Al hops and their activation energies are listed in Table 5. Figure 4 shows the energy profile diagrams for hops A and B. The hop A has the jump distance of 2.97 Å and its activation energy was calculated to be 2.26 eV. A long range diffusion path was constructed using local hops A along the c axis (refer to Figure 3). In this pathway, diffusion of Al 3+ ions is slow. The hop B has a slightly higher activation energy by 0.21 eV than that of hop A. The long-range diffusion pathway connected by hops B exhibit a zig-zag pattern along the c axis.

Self-Diffusion of Aluminium and Oxygen
The ionic conductivity of a material is mainly dependent on the parameters associated with the formation of point defects, pre-exponential factor and activation energy of the migrating ion. The current simulation technique enables the examination of various ion diffusion pathways with their activation energies in Al 2 SiO 5 . Experimental determination of ion diffusion pathways is often difficult. The current methodology has successfully reproduced experimentally determined ion diffusion pathways in a variety of ionic oxide materials [16]. For example, an excellent agreement between the calculated and experimental Li-ion diffusion pathways was observed in LiFePO 4 [40,41].
Vacancy assisted diffusion pathways were examined for Al 3+ and O 2− ions as their Frenkel energies are lower than the other Frenkels. Two different local Al hops (A and B) (refer to Figure 3) with the jump distances of 2.97 Å and 3.14 Å were identified (refer to Figure 3). Individual Al hops and their activation energies are listed in Table 5. Figure 4 shows the energy profile diagrams for hops A and B. The hop A has the jump distance of 2.97 Å and its activation energy was calculated to be 2.26 eV. A long range diffusion path was constructed using local hops A along the c axis (refer to Figure 3). In this pathway, diffusion of Al 3+ ions is slow. The hop B has a slightly higher activation energy by 0.21 eV than that of hop A. The long-range diffusion pathway connected by hops B exhibit a zig-zag pattern along the c axis.      Figure 3.  For the O vacancy diffusion, we identified four local hops. The hop with lowest activation energy (2.75 eV) was considered for constructing a long range diffusion pathway (refer to Figure 5). The oxygen ion involves in a zig-zag type motion along the c axis. A high activation energy of 2.75 eV reveals that the oxygen diffusivity in this materials is low. The other hops exhibited activation energies greater than 4 eV and they are reported here. For the O vacancy diffusion, we identified four local hops. The hop with lowest activation energy (2.75 eV) was considered for constructing a long range diffusion pathway (refer to Figure 5). The oxygen ion involves in a zig-zag type motion along the c axis. A high activation energy of 2.75 eV reveals that the oxygen diffusivity in this materials is low. The other hops exhibited activation energies greater than 4 eV and they are reported here.

Solution of Isovalent Dopants
We considered a variety of trivalent dopants (M = Ga, Fe, Ni, Mn, Sc, Y and La) on the Al site and tetravalent dopants (M = Ge, Ti, Sn, Zr and Ce) on the Si site. Solution energies were calculated using appropriate lattice energies. Buckingham potentials used for the dopants are provided in the supplementary information (refer to Table S1). Isovalent dopants do not introduce any charge compensating defects. Dopants exhibiting low solution energies are useful as those dopants can influence the properties of the material.

Trivalent Dopants
First trivalent dopants were considered on the Al site. The following reaction was used to calculate the solution energy.
Solution energies are reported in Figure 6a. Exoergic solution energy (−0.11 eV) is calculated for Ga indicating that this promising dopant is worth examining experimentally. The possible composition after doping would be (Al1−xGax)2SiO5. The preference of Ga on the Al site is due to the ionic radius of Al 3+ (0.39

Solution of Isovalent Dopants
We considered a variety of trivalent dopants (M = Ga, Fe, Ni, Mn, Sc, Y and La) on the Al site and tetravalent dopants (M = Ge, Ti, Sn, Zr and Ce) on the Si site. Solution energies were calculated using appropriate lattice energies. Buckingham potentials used for the dopants are provided in the supplementary information (refer to Table S1). Isovalent dopants do not introduce any charge compensating defects. Dopants exhibiting low solution energies are useful as those dopants can influence the properties of the material.

Trivalent Dopants
First trivalent dopants were considered on the Al site. The following reaction was used to calculate the solution energy.
Solution energies are reported in Figure 6a. Exoergic solution energy (−0.11 eV) is calculated for Ga indicating that this promising dopant is worth examining experimentally. The possible composition after doping would be (Al 1−x Ga x ) 2 SiO 5 . The preference of Ga on the Al site is due to the ionic radius of Al 3+ (0.39 Å) being close to that of Ga 3+ (0.47 Å). The second promising dopant is Fe 3+ with a solution energy of 0.02 eV. This is supported by the formation of stoichiometric (Al,Fe 3+ ) 2 SiO 5 composition in an experimental study carried out by Grew [11]. Dopants Ni, Mn and Sc exhibit positive (but < 0.45 eV) solution energies, inferring that they can be doped with a small energy penalty. High positive values are observed for Y and La meaning that doping is only possible at high temperatures. This is due to the larger ionic radii of these two dopants as compared to the ionic radius of Al 3+ .

Tetravalent Dopants
Next, the Si site was considered for the doping of tetravalent dopants. The following reaction was used to calculate the solution energy.

MO + Si → M + SiO
Solution energies are reported in Figure 6b. The favorable dopant for this process is Ge 4+ with a solution energy of 0.43 eV. Favorability of the dopant Ge can be attributed to its ionic radius (0.39 Å) being closer to the ionic radius of Si 4+ (0.26 Å) in the tetrahedral coordination. The possible synthesis composition would be Al2Si1−xGexO5. The second most stable dopant is Zr 4+ with the solution energy of 1.17 eV. Solution energies for the other dopants are highly positive meaning that they are unlikely to occur at ambient temperature.

Conclusions
Atomistic simulation based on the classical pair potentials was used to examine the defects, diffusion and solution of isovalent dopants in Al2SiO5 with sillimanite structure. The key defect present in this material is the Al/Si anti-site leading to the formation of disordered configuration as observed in the experiment [11]. The presence of the Al/Si anti-site defect is important as this defect can influence the mechanical and chemical properties of Al2SiO5. The second most stable defect is the O-Frenkel and the concentration of this defect is not significant. Both Al 3+ and O 2− ions exhibit slow

Tetravalent Dopants
Next, the Si site was considered for the doping of tetravalent dopants. The following reaction was used to calculate the solution energy.
Solution energies are reported in Figure 6b. The favorable dopant for this process is Ge 4+ with a solution energy of 0.43 eV. Favorability of the dopant Ge can be attributed to its ionic radius (0.39 Å) being closer to the ionic radius of Si 4+ (0.26 Å) in the tetrahedral coordination. The possible synthesis composition would be Al 2 Si 1−x Ge x O 5 . The second most stable dopant is Zr 4+ with the solution energy of 1.17 eV. Solution energies for the other dopants are highly positive meaning that they are unlikely to occur at ambient temperature.

Conclusions
Atomistic simulation based on the classical pair potentials was used to examine the defects, diffusion and solution of isovalent dopants in Al 2 SiO 5 with sillimanite structure. The key defect present in this material is the Al/Si anti-site leading to the formation of disordered configuration as observed in the experiment [11]. The presence of the Al/Si anti-site defect is important as this defect can influence the mechanical and chemical properties of Al 2 SiO 5 . The second most stable defect is the O-Frenkel and the concentration of this defect is not significant. Both Al 3+ and O 2− ions exhibit slow diffusion with activation energies of 2.26 eV and 2.75 eV, respectively, inferring low ionic conductivity in this material. Exoergic solution energy (−0.11 eV) is calculated for Ga 3+ on the Al site suggesting that the synthesis of (Al x Ga 1−x ) 2 SiO 5 is possible. Ge 4+ is found to be a promising candidate dopant on the Si site though its solution energy is endoergic (0.43 eV). The favorable dopants can be of interest to prevent phase transformation and tune the properties of Al 2 SiO 5 .
Author Contributions: Computation, R.S. and N.K.; writing, N.K and R.S.; analysis and editing, N.K., P.I. and A.C. All authors have read and agreed to the published version of the manuscript.
Funding: This research was financially supported by the European Union's H2020 Program under Grant Agreement no. 824072-HARVESTORE.