The mechanism of Mg diffusion in forsterite and the controls on its anisotropy

Mg diffusion is important for explaining many deformational properties in forsterite but its mechanism is unknown and so the effect of variables such as pressure is difficult to constrain. In this study we used DFT to calculate the diffusion of Mg in forsterite. We predict vacancy diffusion to be highly anisotropic with considerably faster diffusion in the [001] direction while interstitial diffusion is predicted to be more isotropic. Thus we predict that a combination of interstitial and vacancy diffusion is required to reproduce experimentally derived anisotropies. Interstitial diffusion is predicted to be highly pressure dependant such that with increasing pressure the anisotropy of Mg diffusion decreases while temperature has little effect on this anisotropy. Substances like Fe and water likely cause increases in Mg diffusion rate through the creation of extrinsic Mg vacancies and we predict that without modifications to the inherent mobility of Mg vacancies these cause small increases to diffusional anisotropy at 1300 and 1600 K but very large increases at 1000 K.


Introduction 42
Diffusion of cations occupying the octahedral metal sites in olivine controls processes that are active 43 in the Earth's crust and upper mantle, and which underpin a range of geophysical and geochemical 44 techniques. In the upper mantle, where olivine with composition close to (Mg0.9,Fe0.1)2SiO4 is the 45 dominant phase, the diffusivity of Mg is important in understanding electrical conductivity (Fei et al., 46 2018, Yoshino et al., 2009, Yoshino et al., 2017, Schock et al., 1989) and could influence deformation 47 even though Mg is a rapidly diffusing species as argued in Jaoul (1990). Anisotropic Mg diffusion could 48 be an important factor in explaining the anisotropic conduction seen in high conductivity layers 49 underneath young oceanic plates (Fei et al., 2018) and, if Mg diffusion is important in forming olivine 50 textures, could also help explain the variety of textures that are formed by olivine under different 51 conditions (Karato et al., 2008). Mg-Fe interdiffusion occurring in zoned phenocrysts from volcanic 52 products is increasingly used as a petrological tool (diffusion chronometry) to understand the 53 timescales of pre-eruptive processes operating in the days and weeks prior to eruption (e.g. Mg and Fe in olivine and spinel can be used to infer the post-crystallisation thermal history of 56 ultramafic igneous bodies (Ozawa, 1984). Diffusion can also lead to magnesium and iron isotope 57 fractionation (Teng et al., 2011). 58 Our understanding and ability to model all of these processes relies on accurate determination of the 59 Mg self-diffusion and Fe-Mg interdiffusion coefficients in olivine and thus this has been the focus of a 60 range of experimental and computational studies reviewed by Chakraborty (2010). However, details 61 of the atomic scale basis of Mg self and inter-diffusion in olivine have thus far eluded a full atomistic 62 explanation and this limits our ability to confidently make use of this data under the wide range of 63 conditions where diffusion is important. In this work we shall study the atomistic mechanisms of Mg 64 self-diffusion as the more straightforward of these two processes. 65

Experimental measures of Mg diffusion
Previous experimental studies have identified several key features of Mg self-diffusion as well as 67 questions that remain unanswered. Despite early uncertainty, it is clear that magnesium self-diffusion 68 is faster than the self-diffusion of oxygen or silicon (for a review of this history see Chakraborty (2010)). 69 In detail, magnesium self-diffusivity is found to be mildly sensitive to pressure, to be anisotropic and 81 to depend on the chemistry of the olivine crystal. Diffusion along [001] has been found to be faster 82 than diffusion along [100], which is faster than diffusion along [010] (Chakraborty et al., 1994) though 83 other studies have found diffusion along [010] to be faster than diffusion along [100] (Andersson, 84 1987, Jollands et al., 2020). There are some differences in experimental activation volumes 1-3.5 85 cm 3 /mol (Chakraborty et al., 1994) or 4.0-4.6 cm 3 /mol (Fei et al., 2018) but in all cases these are small 86 and so pressure has little effect on diffusion rates. These activation volumes come from data solely in 87 the [001] direction which is important because they do not reflect all processes that occur in the crystal 88 as shall be explored in the text. 89 Although these experiments provide the critical data needed to model diffusion-controlled processes 90 in olivine, several aspects of Mg diffusion remain enigmatic and some parameters have not been fully 91 established. For example, the reason for the anisotropy of diffusion is not clear and the effect of vacancies which are most stable on M1 rather than M2 sites and octahedrally coordinated Mg 119 interstitials which form a split-interstitial structure (two magnesium ions in tetrahedral coordination 120 located on opposite sides of the M1 site) is stable (Walker et al., 2009). The mobility of some of these 121 defects has been studied using interatomic potentials (Bejina et al., 2009, Jaoul et al., 1995, Walker et 122 al., 2009 where it was found that Mg vacancies are more mobile than Mg interstitials (Walker et al., 123 2009), that pressure has a limited effect on mobility along the M1 chain as was found in experiments 124 (Jaoul et al., 1995, Bejina et al., 2009  Compared to the timescale accessible to direct atomic scale simulation using molecular dynamics, 137 point defect diffusion in minerals is usually slow. Methods available to simulate diffusion thus seek to 138 describe diffusion by repeated rare events which can be studied in detail, and then combined in order 139 to describe diffusion on a meaningful timescale. The rare events are typically hops of point defects 140 between adjacent sites. For example, one of a number of atoms could migrate into a vacancy, 141 effectively moving the vacancy and permitting diffusion via a vacancy mechanism, or an interstitial 142 atom could move into one of a number of different interstitial sites, permitting diffusion via an 143 interstitial mechanism. Repeated occurrences of these hops leads to a random walk of the defect and bulk self-diffusion (Tilley, 1987). Our approach to simulating Mg diffusion in forsterite thus follows 145 three steps. First, we make use of density functional theory to determine the structure and relative 146 stability of stable Mg point defects in forsterite. These models represent the ground state end-points 147 of the hops leading to diffusion. Second, we probe the energy landscape that must by traversed by 148 the defect during a hop. This provides us with the energy barrier that must be overcome for the hop 149 to proceed and the structure of the transition state (the configuration with maximum energy on the 150 minimum energy pathway between the start and the end point We created models of Mg vacancies by removing an Mg 2+ ion from an M1 or M2 site in a 171 (2x1x2) forsterite super cell. Interstitial defects were created by inserting an extra Mg 2+ ion into 172 potential interstitial sites in the structure. In both cases the cell parameters were always fixed to those 173 of the defect free crystal to approximate the dilute limit. To account for atomic relaxation around the 174 defects, the structure was then relaxed until the forces on all atoms were less than 0.01 eV/Å and an 175 energy change between different geometric steps was less than 110 -5 eV/atom. Repeating 176 calculations with increased cutoffs changed the energy of the supercell by <0.1 meV/atom. A (2x1x2) 177 forsterite supercell was used to ensure that there was roughly 10 Å between repeating vacancies in 178 all directions, a distance we found to be sufficient to contain the important atomic relaxations. 179 Simulation cells containing vacancies or interstitials have a net charge and so the energy calculated by 180 CASTEP includes a defect-defect interaction term between adjacent supercells which does not reflect 181 our desired energy of a charged defect in an infinite medium. We can approximately correct for this 182 interaction by assuming it is the energy of a periodic array of point charges in a uniform neutralising 183 background charge. This was done using the method of Leslie and Gillan (1985), first used for 184 forsterite by Brodholt (1997). To use this method the relative permittivity of the cell needs to be set 185 -we used a value of 6.2 (Weast, 1981). We repeated these calculations for a (4x2x4) supercell 186 containing a Mg vacancy, and the vacancy energy changed by <0.01 eV, suggesting that our simulation 187 cell size and energy corrections are sufficient for our needs. Once ground state structures and energies for the defects had been determined, we enumerated the 194 possible hops (where a defect moves from location to location) and for each hop we determined the 195 pathway and found the transition state structure and energy. We did this by using a constrained optimisation approach. We first determined an approximate path for the hop (for vacancy diffusion 197 this consists of two vacancies with a Mg atom at a point between the vacancies, for interstitial 198 diffusion the interstitial atom is located between stable interstitial sites). For each hop we tried 199 multiple paths, but direct paths proved to have the lowest transition state energy in all cases. state. While this method may not definitely find the transition state our frequency calculations 207 (below) typically returned a single imaginary eigenvalue of the dynamical matrix, as expected for a 208 transition state. In the few cases, which were all for interstitial diffusion, where this was not the case 209 the candidate transition state was found by manual adjustment based on visualising the eigenvectors 210 of the imaginary phonon frequencies until a single imaginary eigenvalue was found. It turned out that 211 this manual adjustment changed the activation energy of the hop by <0.01 eV suggesting that the 212 constrained optimisation method is highly reliable for finding activation energies even if they are in 213 complex parts of the energy hypersurface. 214 We repeated the calculations described above at 0, 5 and 10 GPa and 1000, 1300 and 1600 K by setting 215 the simulation cell dimensions to minimise the Gibbs free energy of the defect free cell. The effect of 216 pressure is easily accounted for by adding the PV term to the internal energy of the system. The effect 217 of temperature requires consideration of the thermal motion of the atoms. We include this effect by 218 making use of lattice dynamics to evaluate the phonon frequencies and then use these to evaluate 219 the vibrational entropy of the crystal. Phonon frequencies were determined using the finite 220 displacement method of CASTEP with finite displacements of 0.01 bohr. All lattice dynamics 221 calculations were performed solely at the q=(0,0,0) point. While this calculation at a single q-point may introduce a significant sampling error all of our calculations involve comparisons between two 223 very similar structures -the start/end point of a diffusion step and its transition state -and so the 224 effect of sampling errors are likely to be small but this is a limitation of the method. For lattice 225 dynamics we tightened the convergence criteria on the forces and energy for the geometry 226 optimisation to 0.001 eV/Å and 110 -9 eV/atom, respectively. A few end points and transition states 227 were sampled with 0.00075 eV/Å and 5x10-10 eV/atom cuts off and the change in free energy caused 228 by these increased cutoffs was <1 meV/atom. We determined the Gibbs free energy at a wide range 229 of temperatures and at least 5 different volumes and then the energy at each volume with the 230 following equations: 231

Equation 3 234
Where U(V) is the internal energy and νk,i(V) is the frequency of the phonon with wave vector k in the 235 i-th band at volume V. At the pressure and temperature of interest the appropriate volume and 236 energy was determined by fitting 2 nd order polynomials across our volume range and minimising 237 Equation 1. This method is quasi-harmonic as it ignores anharmonic effects beyond those caused by 238 thermal expansion. 239

From defects to diffusion 240
The self-diffusion of a Mg by a vacancy mechanism can be represented by: 241

Equation 4 242
Where is the diffusion coefficient of Mg vacancies and NVac is the atomic fraction of Mg 243 vacancies. 244 As shown below, our atomic scale simulations suggest that diffusion of both interstitials and 245 vacancies can be important for magnesium diffusion in pure forsterite. To account for this possibility we use the assumption that vacancies and interstitials diffuse independently of each other, which 247 means that the total self-diffusion of Mg in forsterite is given by: 248 Other diffusing species (which are not considered in this paper) would have their own term if present. 250 For systems with simple geometry, the diffusion coefficients can be found analytically from the 251 attempt frequency, the migration entropy, the activation energy and the crystal structure. For 252 example, for a single hop the coefficient is given by (Poirier, 1985): 253 where α is a geometric prefactor to account for the degeneracy of the hop, q is a dimensionality 255 constant (q = 2, 4 or 6 for 1, 2 or 3D diffusion), l is the length of the hop and the two exponential terms 256 are the migration entropy and the migration enthalpy, respectively. This approach has been used to For our KMC method we need to know the concentration of defects and the rate at which each Mg 264 hop can occur. As explained above the concentration of intrinsic defects (vacancies and interstitials) 265 were determined by minimising the free energy of the Frenkel reaction at the appropriate P and T. To 266 determine the rate of hopping we used lattice dynamics to probe the vibration of atoms around the 267 point defects in their ground state and transition state configurations. This allows us to model the 268 effect of temperature on point defect mobility. The rate, k, at which a defect hops from one location 269 to another is given by: 270 where v is the attempt frequency (in Hz). The activation energy term was calculated from our 272 constrained optimisation. In order to calculate the attempt frequency and activation entropy we used 273 Vineyard theory (Vineyard, 1957) which is based on absolute rate theory. Both of the temperature-274 based factors (vibrational entropy and attempt frequency) are combined into a modified attempt 275 frequency (v*) which is found from the ratio of the calculated phonon frequencies: Where q is the dimensionality constant as above. 313 3 Results 314

Defect Energies and Concentrations 315
There are two sites for Mg vacancies in forsterite -the M1 and the M2 sites. We calculate that M1 316 sites are strongly favoured over M2 sites with pressure increasing the vacancy preference for M1 sites 317 (Table S1). This preference for M1 over M2 vacancies agrees with previous calculations though there 318 is some difference in the energy of this preference (0.9-1.2 eV in this work, ~1.9 eV with forcefield 319 calculations (Walker et al., 2009) or ~0.8 eV previously using DFT (Brodholt, 1997)).
We have also considered Mg interstitials. As with Walker et al. (2009) we found that the most stable 321 position is a split interstitial at the M1 site with 2 Mg atoms displaced from the centre of this site in 322 opposite [010] directions (shown in Figure S1). This arrangement is very stable with alternative 323 arrangements of the Mg at this site all relaxing into this one. Even placing a Mg atom in an I1 site 324 causes it to relax into this split interstitial arrangement. The other stable configuration is found by 325 placing an additional Mg in the I2 site. The Mg interstitial in the I2 site has an octahedral coordination 326 like the M1 and M2 and is thus geometrically similar to them. At 0 GPa the split M1 interstitial is 327 slightly favoured over the I2 arrangement (~0.2 eV) but with increasing pressure the I2 configuration 328 is favoured (Table S1) as the split M1 arrangement is larger than the I2 arrangement. In QM-MM 329 embedded cluster calculations (Walker et al., 2009) the split M1 geometry was found to be favoured 330 over an I1 interstitial geometry by ~4.4 eV but an I2 geometry was not reported. In our own forcefield 331 calculations we were unable to stabilise an I2 arrangement as I2 arrangements always relaxed into M1 332 arrangements. Forcefields are thus likely poor at calculating these interstitial structures. where a is a reaction vector for the Frenkel reaction (between 0 and 1), ΔE is the energy of the Frenkel 346 reaction and Sconfa is the configurational entropy after the reaction has proceeded forward by a. The 347 results of this minimisation are given in Table 1 To simulate the other defects we probed all likely sites, found the sites with the minimum enthalpy 360 and then calculated their high temperature energy through Equation 1-3. The energy of these 361 reactions is shown in Table S2 but all reactions have substantially higher energies than the Frenkel 362 reaction (R1). By including these other reactions, which also compete in configurational entropy 363 space, we change the concentration of Mg defects by less than 0.001% and thus these can safely by 364 ignored and we shall only consider the Frenkel reaction R1 from now on. 365 366

Vacancy Hops 367
For Mg diffusion by vacancy hopping we found six different vacancy diffusion hops for which we 368 calculated the geometries and energies of hopping. The hops that we have considered are shown and 369 labelled in Figure 1 with their dimensions listed in Table S3 and described in the SI.
The activation energies and frequencies of these hops are presented in Table 2 and the barriers to  371 diffusion are shown in Figure 2. Notably the A hop which is directly along the [001] direction has a 372 substantially lower activation energy than all other M1 hops. The easiest hop from an M2 site is the 373 C hop back to an M1 site. These two effects combine such that vacancies will diffuse easily along the 374 [001] direction when in a M1 site and will have difficulty escaping to an M2 site. If they do escape to 375 an M2 site they will be converted quickly back to an M1 site. The weighted probability of these hops 376 is shown in Figure 2 and an alternative representative in Figure S2  was an effect of simply using DFT as against using forcefields we recalculated our results using GULP 382 with a TBH1 forcefield (Wright and Catlow, 1994) (Table S4, computational details in supplementary 383 information). We find that generally DFT produces lower barriers than forcefield calculations but that 384 the order of the hops is the same with both DFT and forcefield calculations. Crucially the activation 385 energy of the easiest A hop (which largely controls the overall diffusion) is very similar with both 386 methods 0.77/0.75 eV which means that both DFT and forcefield calculations return a very similar 387 diffusivity for anhydrous vacancy diffusion. 388 We also considered the effect of pressure on the activation energies of these vacancies. As shown in 389 Table S5 and Table S6 going from 0 to 10 GPa makes negligible differences to the activation energy or 390 * of any of the hops. The small differences seen are miniscule compared to the effect pressure has 391 on the vacancy concentration as described above. 392

Mg interstitial hops 393
As Mg interstitials occupy M1 and I2 sites-the latter of which are simply shifted M2 sites-the relative 394 geometry of interstitial hops are identical to those of vacancies. These hops are pictured and labelled 395 in Figure 3 and their barriers in Figure 4 (and tabulated in Table S7) with their energies and frequencies in Table 2 (and more pressure derivatives listed in Table S6 and S8). The probability of any of the hops 397 occurring is shown in Figure 4 and alternatively in Figure S3. These are again described in the 398 supplementary information. 399 Interstitial hops I and J, which are between M1 and I2 sites, are the most favourable with activation 400 energies <0.6 eV. In part this is because in the split M1 configuration one Mg at the M1 site is already 401 close to an I2 site. Pressure has a small effect on the attempt frequency (Table S6) but a relatively 402 large effect on the activation energy of these hops (Table S8)

Diffusion 408
Using our KMC algorithm we can convert hops into diffusion rates. The diffusion coefficients for both 409 vacancy and interstitial hopping are presented in Table 3 (these are listed at 5 and 10 GPa in Table S9  410 and S10). Vacancy diffusion is highly anisotropic with diffusion along [001] being orders of magnitude 411 faster than diffusion along [100] or [010]. This is an outcome of diffusion where the hop directly along 412 [001] is ~0.75 eV easier than any other M1 hop. In the absence of any additional undiscovered 413 hops/mechanisms this will always hold. Interstitial diffusion is much more isotropic than vacancy 414 diffusion due to the favourability of M1 to I2 hops (I and J) which go in all three primary directions. 415 To calculate total diffusion of Mg in forsterite we added together the rates of Mg vacancy and 416 interstitial diffusion. This assumes that Mg Frenkel pairs are not associated with each other. To test 417 this assumption, we calculated the binding energy of this pair by running separate simulations with 418 isolated Mg vacancies and interstitials and then calculations with them adjacent in the same unit cell 419 and comparing the difference in enthalpy. We find that the binding energy is approximately -1.9 eV 420 with a negative number indicating that bound defects are more stable than unbound defects. This is 421 a large number but it is much smaller than the configurational energy gains of randomly scattering Mg vacancy and interstitial pairs for low concentrations. For the pairing energy to exceed this 423 configuration entropy, the defect concentration would need to be above 1.2x10 -3 defects per unit cell 424 at 1300 K, many orders of magnitude larger than the predicted vacancy concentrations (Table 1) Table S11) as a function of pressure. Notably we find a larger pressure derivative for 458 intrinsic diffusion coefficients than has been seen in the literature (Chakraborty et al., 1994, Fei et al., 459 2018. Our activation volumes are 6.69 cm 3 /mol at 1000 K, 7.51 cm 3 /mol at 1300 K and 7.84 cm 3 /mol 460 at 1600 K. The pressure dependence of diffusion is strongly controlled in our calculations by the 461 pressure dependence of defect concentration (Table 1) with little effect of the defect mobility (Table  462 3

Anistropic intrinsic diffusion 475
One of the most notable features of our results is that Mg diffusion can be strongly anisotropic. Figure  476 7 shows the anisotropy of this diffusion as a function of pressure. We find that anisotropy decreases 477 with pressure due to the increasing importance of interstitial diffusion, which is less anisotropic, while 478 temperature has little effect on anisotropy. At 1600 K and 0 GPa (corrected) we find the ratio of

Anisotropy Changes in the Upper Mantle 502
While the dependence of anisotropy on pressure is large this probably has little implication in the 503 upper mantle. After applying pressure corrections a 0-10 GPa range in the upper mantle would be 504 equivalent to ~4-16 GPa in our pressure scales. The largest changes in anisotropy come at the lowest 505 pressures and so across the pressure range of the upper mantle, changes in Mg diffusional anisotropy 506 with depth will typically be up to an order of magnitude except at the coldest temperatures (1000 K) 507 where this could reach 1.5 orders of magnitude. These changes are likely to too small to have any 508 major effects on mantle rheology that change with depth. 509  (Chakraborty, 2010). In the Fe-Mg interdiffusion case the mobility of vacancies is also 522 modified as they include Fe self-diffusion coefficients. In the case of water the mobility of (2 ) 523 could be different to ′′ but likely similar. In these and other cases we expect the change in the 524 concentrations of vacancies to generally outweigh the changes to the mobility of vacancies due to the 525 small number of intrinsic defects produced by the Frenkel reaction (Table 1). Extrinsic Mg vacancy 526 concentrations can be many orders of magnitude higher than our predicted intrinsic Mg vacancy 527 concentration in many systems. This is seen in iron-containing olivine (Dohmen and Chakraborty, 528 2007) where the Mg vacancy concentration is many orders higher than predicted here due to R7. Thus 529 the prime reason that various contaminants cause an increase in Mg diffusion rates is likely to be the 530 production of more Mg vacancies. Critically Mg vacancies can be produced in this way but producing 531 extrinsic Mg interstitials is much more difficult. This means that extrinsic defects are likely to produce 532 a strong imbalance in the Mg vacancy vs Mg interstitial ratio. (1-100 ppm) (Fei et al. 2018). At 1300 K we predict these defects to increase diffusional anisotropy 547 (compared to perfect forsterite) by 2-5 times at 5-10 GPa (uncorrected). As temperature increases 548 this effect decreases such that by 1600 K iron and water increase diffusional anisotropy by less than predicted rates at 0 GPa (corrected-see supplementary information) determined by fitting between 706 our pressure corrected values (the same plot with a 5 GPa pressure correction is shown in Fig S4)  direction) for different diffusion rates (Dx) in a system of self diffusion+extrinsic vacancies. Anistropy 826 is relative to pure forsterite which is 1. This was determined by solving Equation 13 as a function of 827 diffusion rate. Lines are at 1000 (blue), 1300 (orange) and 1600 K and solid lines represent 5 GPa, 828 dotted lines 10 GPa (uncorrected) which correct to around 1 and 6 GPa respectively. 829 The dark region represent the range of Dx between Fe=1-20% for olivine at 0 GPa, 1300 K and 830 fO 2 =10 -7 Dohmen and Chakraborty (2017). The ligt region represents the range of Dx for water 831 ranging between 1-150 wt. ppm at ~8 GPA and 1300 K (Fei et al. 2018 Table 2: Activation energy and modified attempt frequency v* of various hops (shown in Figure 1  846 and 3 with the hop distances outlined in Table S4 and S8) for hydrous and anhydrous forsterite at 0 847 GPa uncorrected. Hop L could not be stabilised but is very high in energy. Hops with an asterisk go in 848 the reverse direction where this is not equivalent. 849