Atomistic simulations of Mg vacancy segregation to dislocation cores in forsterite

Interactions between dislocations in olivine and extrinsic cation vacancies created under hydrous or oxidizing conditions may influence the rheology of the Earth's upper mantle. In this study, we use atomic-scale simulations to calculate segregation energies for bare and protonated Mg vacancies to M1 and M2 sites in the core regions of [100](010) and [001](010) edge dislocations, and [100] and [001] screw dislocations. Calculated segregation energies are different for the two symmetry distinct M sites. The segregation energies calculated for the tightest binding M1 sites around [100] screw and [100] (010) edge dislocations are comparable to those calculated for the tightest binding M2 sites. Concentrations of M2 vacancy-related defects will thus be low in the core regions of these dislocations, given the comparatively high energy of these defects in the bulk lattice. In contrast, segregation energies for M2 defects to [001](010) edge dislocation cores are considerably lower than for equivalent M1 defects, and M2 vacancy concentrations around these dislocations will be similar to M1 vacancy concentrations. This means that the effect of magnesium vacancies on the mobility of the [001](010)


Introduction
Although nominally anhydrous, under the pressure and temperature conditions of the Earth's mantle, olivine can incorporate modest quantities of water into its crystal structure, primarily as protonated cation vacancies (e.g. Martin and Donnay, 1972;Bai and Kohlstedt, 1993;Kohlstedt et al., 1996). The water solubility limit in olivine is sensitive to the water and oxygen fugacities, and also increases with silica activity, consistent with incorporation via protonation of M site vacancies (Gaetani et al., 2014).
Protonated vacancies interact with dislocations in the olivine crystal lattice, and may preferentially occupy atomic sites near such defects over sites in the bulk lattice, with the effect that vacancy-related defect concentrations are potentially much greater in the core region of a dislocation than in the unstrained bulk lattice. In extreme cases, impurity segregation can affect the bulk chemistry of minerals, as for example in the formation of striped chemical zoning in olivine during low strain-rate deformation, attributed to Fe 2+ segregation to sub-grain boundaries formed by aligned arrays of edge dislocations (Ando et al., 2001). High-resolution synchrotron images show that the concentration of protonated defects in olivine is greatest around grain boundaries and cracks, demonstrating that these defects in olivine will tend to segregate to locally strained regions of the crystal lattice (Sommer et al., 2008), while hydrogen concentrations near [001] dislocations in water-saturated olivine can be sufficiently great to induce climb dissociation of the dislocation core (Drury, 1991). cores, pinning the dislocation in place and reducing strain rates accordingly (Cottrell and Bilby, 1949).
However, some defects can increase dislocation mobility in the glide creep regime by reducing the Peierls barrier to glide, including vacancies in fcc Al (Lauzier et al., 1989;Lu and Kaxiras, 2002), hydrogen in Fe (Taketomi et al., 2008), or interstitial O defects in UO 2 (Ashbee and Yust, 1982;Keller et al., 1988). In olivine, interactions between water-related defects, occurring primarily as protonated vacancies may facilitate deformation by increasing dislocation mobility in the dislocation climbcontrolled creep regime (e.g. Mackwell et al., 1985;Chen et al., 1998;Girard et al., 2013). Hydrous defects are thought to reduce the Peierls stress, s p , required for dislocation glide, whose measured value for hydrated olivine is ~1.6-2.9 GPa (Katayama and Karato, 2008), considerably lower than values measured for dry olivine polycrystals, which range from at least 3.8 GPa (Idrissi et al., 2016) to as much as ~15 GPa (Demouchy et al., 2013). Recent forced-oscillation measurements have suggested that Mg vacancies, produced to charge balance the oxidation of Fe 2+ to Fe 3+ (Stocker and Smyth, 1978;Nakamura and Schmalzried, 1983) may enhance attenuation in Fo 90 olivine (Cline et al., 2018).. Natural dunites deforming in the dislocation creep regime also show a moderate sensitivity to the oxygen fugacity (Keefner et al., 2011). This suggests that bare cation vacancies may have a similar influence on the mechanical properties of olivine as protonated cation vacancies.
The short length scales characteristic of impurity segregation to dislocation cores mean that it can be difficult to study experimentally, although developments in the field of atom probe tomography mean that it is now possible to visualize impurity clouds around dislocation lines (Miller et al., 2006;Peterman et al., 2016). Theoretical modeling offers an alternative approach, allowing direct access to the atomic scale and control over system chemistry. While interactions between dislocations and point defects far from the dislocation line can be adequately modeled using linear elasticity theory, is the dislocation core non-elastic, atomic-scale relaxation can be substantial. One way to model a dislocation is to insert two or more dislocations into a 3D-periodic simulation cell, with their Burgers vectors b summing to zero to ensure continuity at the boundaries. Although this cell can be sufficiently small as to make the use of ab initio methods practical, care must be taken to minimize dislocation-dislocation interactions. An alternative is to embed a single dislocation in an isolated cluster of atoms with periodic boundary conditions along the dislocation line (Walker et al., 2005a). Both the cluster-based (Walker et al., 2005b) and supercell (Mahendran et al., 2017) approaches have been used to simulate [100] and [001] screw dislocations in forsterite, producing comparable dislocation core structures, although the latter study did not report core energies for either dislocation.
While atomistic modeling is a powerful tool for studying dislocations and their interactions with point defects, there are several limitations that restrict its range of applicability. Firstly, obtaining a converged dislocation core structure and energy may require the use of very large simulation cell, containing many hundreds or thousands of atoms, for which the computational cost of using quantum chemical methods such as DFT can be prohibitive. Instead, as in this study, interatomic potentials are more commonly used, which are parameterized by fitting to experimental data or ab initio calculations. The second problem is that the dislocation itself breaks the translational symmetry of the crystal, meaning that interactions between point defects cannot be parametrized using any of the techniques available for solid solutions, such as cluster expansion (Sanchez et al., 1984) or Special Quasirandom Structures (Zunger et al., 1990), and the dislocation energy must be obtained from fully atomistic calculations. In practice, this limits calculations to the dilute limit.
In this study, we use the cluster-based approach to determine segregation energies for bare and protonated cation vacancies to dislocation cores in the forsterite (Mg 2 SiO 4 ) end-member of the olivine solid solution. Since the silica activity in mantle peridodites is high, Mg vacancies are expected to be more abundant than Si vacancies. Consequently, we consider only Mg vacancies on the two symmetry distinct M sites, labeled M1 and M2. At low pressure, the easiest slip system in olivine is [100](010), but the [001](010) slip system becomes more active at high pressure (Couvy et al., 2004;Hilairet et al., 2012 (Kröger and Vink, 1956), to sites within the core regions of these dislocations. Energies for the different defects are compared to elucidate the effect of site occupation and hydrogen fugacity on the interaction between Mg vacancies to dislocations. In the bulk, M1 vacancies are energetically more favorable than M2 vacancies in the bulk lattice (Brodholt, 1997), but this may be different near dislocation cores, which may have implications for olivine rheology in dislocation-controlled creep regimes, as well as for Mg diffusion in crystals with high dislocation density.

Computational methods
Dislocation core structures and segregation energies were calculated using the cluster-based approach, in which an isolated dislocation is inserted at the axis of a 1D-periodic cylinder of atoms (Sinclair, 1971;Walker et al., 2005ab). The starting coordinates for the atoms are determined from the elastic displacement field u(r) calculated using the sextic formulation for a dislocation in an anisotropic medium (Stroh, 1958). For edge dislocations, this is a non-conservative algorithm and atoms must be removed from the simulation cell to obtain a physically reasonable initial dislocation structure. To do this, a branch cut is created that is normal to both the Burgers and dislocation line vectors. Any atoms that are displaced across this branch cut by the displacement field u(r) are deleted. Atoms in close proximity to the branch cut are merged with any nearby atoms, if the distance between them falls below a specified threshold d min . The cluster of atoms is subsequently divided into two concentric regions, with radii R I and R II . During the geometry optimization step, atoms in the inner region (region I) are permitted to relax freely, while those in the outer region (region II) are held fixed at the coordinates predicted using the elastic displacement field.
The total excess energy per unit length, E dis , contained within radius r of an isolated dislocation is where K is the elastic energy coefficient and depends on the dislocation geometry and elastic constants C ij , E core is the energy contained within the core region (termed the core energy), and r c is the radius of the dislocation core, within which the displacement field diverges from the predictions of linear elasticity. The core radius r c is an undetermined parameter, whose value cannot determined from the radial excess energy of the dislocation. Its value must be chosen in order to set a gauge for the core energy. In this study, we use a core radius of 2b, where b is the absolute magnitude of the Burgers vector.
The core energy is determined from atomistic cluster-based simulations by fitting equation (1) to the calculated radial dependence of the excess energy, which is the difference between the energy of a cluster containing the dislocation and a reference system containing an identical number of atoms. E core is also the excess energy of the dislocation at r = r c . The excess energy is calculated from the energies of the individual atoms as where E dis (r) is the total energy of the atoms within r of the dislocation line, the sum runs over the different atomic species present, n species (r) gives the number of atoms of each species within r, and E species is the energy of the species in the bulk lattice. This is equal to where E vac is the energy of a supercell from which one atom of the specified type has been removed, without relaxing the coordinates of the remaining atoms, E supercell is the energy of the supercell without a vacancy, and E isolated is the energy of an isolated atom of the specified type. In single-component crystals, this is identical to the energy of the unit cell divided by the number of atoms it contains.
The core energy and core displacement field of a dislocation in a two-region cluster depends on the radius of the relaxed region. A region I radius R I = 25 Å was sufficient to guarantee convergence of the which uses a charge neutralizing term to guarantee convergence of the energy at a finite distance. A cutoff range of r cut = 15 Å and damping parameter x = 0.2 Å -1 were used, giving lattice parameters and elastic constants that differ from the values calculated using the Ewald method by <1%. As electrostatic interaction between ions is truncated at r cut , the region II radius R II of R I +r cut is used for all cluster calculations.
Due to the large size of the simulation cell, all calculations are performed using empirical interatomic potentials in the program GULP (Gale, 1997;Gale and Rohl, 2003). The interatomic potentials used are from the THB1 model, which was parameterized by fitting to experimental data (Sanders et al., 1984; Lewis and Catlow, 1985), and reproduces the physical properties of forsterite with reasonable accuracy (Price et al., 1987). Following Wright and Catlow (1994), we model protonated vacancies using the parameters developed by Schröder et al. (1992) to treat (OH)groups in zeolite, incorporating the subsequent modifications made to the Morse potential by Gatzemeier and Wright (2006). This potential, labeled THB1, has been widely used to model point and extended defects in forsterite, including Mg point defects (Walker et al., 2009), surface structures and energetics (de Leeuw et al., 2000), and screw dislocation core structures and energies (Walker et al., 2005b).
In cluster calculations, the segregation energy E seg of a single point defect at an atomic site in a dislocation core is determined by calculating the excess energy DE dis of a point defect of the specified type embedded it in a simulation cell whose length is a multiple of the unit cell edge parallel to the dislocation line vector, x, and comparing it with the excess energy DE perf of an isolated defect in the bulk lattice, taken here to be the excess energy DE perf of a point defect embedded in a 3D-periodic supercell of the material. This is equivalent to where E dis is the energy of a cluster containing a dislocation, E dfct+dis is the energy of that same cluster with a single point defect inserted, E supercell is the energy of a defect-free 3D-periodic supercell, and E dfct+supercell is the energy of a supercell containing a point defect. Negative segregation energies indicates that the point defect will tend to bind to this site to lower the total energy of the system, while positive segregation energies indicate the reverse.

Dislocation core properties
The energy of a given dislocation depends on its coordinates within the crystallographic plane normal to x. For each of the dislocations considered in this study, there are several possible symmetrically distinct origins (labeled in Fig. 1). In the case of edge dislocations, which also break any rotational symmetry of the crystal about the line vector x, the number of symmetrically distinct origins for a dislocation can be even higher. Core energies for the most stable configuration found for each (1) are reported in Table 1. Also shown are their associated elastic energy coefficients K, which are determined from the elastic constant C ij using the Stroh sextic theory (Stroh, 1958 However, the core radius r c depends on the length of the Burgers vector, which is shorter for the [100] (010) edge dislocation, and E core therefore corresponds to the energy of a smaller region. As can be seen in  [001](010) edge dislocation is symmetric, due to the existence of mirror planes parallel to (001) located at z = 0.25 and z = 0.75, passing through the row of Si atoms parallel to [010]. Both edge dislocations lie on the median planes of the M2O 6 polyhedra (y = 0.25/0.75). This is consistent with quantum mechanical calculations of generalized stacking fault energies, which find that ideal shear stresses for [100](010) and [001](010) slip are lowest when slip is localized at y = 0.25 (Durinck et al., 2005). The Peierls stresses for dislocations gliding on (010) are similarly lowest when glide is on the plane at y = 0.25 (Durinck et al., 2007).
For the [100] screw dislocation, we find that the origin of the most stable core structure is (0.5, 0.25), halfway between adjacent M2 sites (labeled site C in Fig. 1a), which has a calculated core energy of E core = 1.78 eV/Å. This core structure was also reported by Mahendran et al. (2017), who used the alternative supercell approach. Earlier work using the cluster-based approach, by contrast, found that the dislocation centered on the M1 site has a lower energy (Walker et al., 2005b), for which we compute a relatively high core energy of 1.97 eV/Å. The discrepancy is likely due to the fact that Walker et al. (2005b) searched for the minimum energy core structure using single point energy calculations at each possible core position, whereas the core structures were relaxed in this study. Local atomic-scale structure thus has a determining effect on the relative stability of the different core configurations for the [100] screw dislocation in forsterite.
Whereas other dislocations gliding on (010) are located on the median plane of the sheet of M2O 6 octahedra, the most stable core structure of the [001] screw dislocation is centered on the column of M1O 6 polyhedra running parallel to [001] (labeled site D in Fig. 1b), consistent with previous theoretical calculations (Walker et al., 2005b;Mahendran et al., 2017). As found in previous studies (Carrez et al., 2008), the [001] screw dislocation has a non-planar core. This is can be seen clearly in where z is the coordinate along the dislocation line, to disappear. In what follows, the region in which 0.0 <= z < 0.5 is referred to as the "lower" region, and the region with z satisfying 0.5 <= z < 1.0 as the "upper" region. In this labeling scheme, the lower region corresponds to those M1 sites that relax away from the (010) glide plane, and the upper region to the sites that relax toward it.

Excess energies of defects in the bulk lattice
Segregation energies are calculated from equation (4) Creating an M1 vacancy, whether protonated or bare, is thus more energetically favorable than creating an M2 vacancy. DH M1→M2 , the enthalpy required to exchange an Mg vacancy between the M1 and M2 sub-lattices is 1.9 eV, identical to previous values of DH M1→M2 calculated using empirical potentials (Jaoul et al., 1995;Walker et al., 2009) calculations (Brodholt, 1997). The energy difference DH M1→M2 between the {2H M1 } X and {2H M2 } X defects, at 2.4 eV, is even greater than that for bare vacancies. As the relative concentrations of vacancies on the two sites depends exponentially on DH M1→M2 , M1 vacancies, whether bare or protonated, will be considerably more abundant than similar M2 vacancy-related defects in the bulk lattice.

Segregation of M1 vacancies
The segregation energy for the {V M1 }″ defect around a [100](010) edge dislocation is lowest for the three sites located directly below the dislocation line (Fig. 4a). The {V M1 }″ defect binds particularly tightly to the site directly below the glide plane, which has a segregation energy of -3.00 eV. which is located directly below the glide plane and on either side of the dislocation line, E seg = -1.08 eV, and segregation energies are only marginally higher for M1 sites above the glide plane.
The calculated minimum segregation energies for M1 vacancies binding to screw dislocation cores are higher than those for the edge dislocations, consistent with the lower stresses induced by a screw dislocation. For the [100] screw dislocation, the low energy sites are distributed radially around the dislocation core (Fig. 6), with the tightest binding sites being those closest to the (010)

Segregation of M2 vacancies
The lowest segregation energy site for {V M2 }″ around the [100](010) edge dislocation is not at the glide plane, but at x = 0 on the first sheet of M2O 6 octahedra below the dislocation (Fig. 4) (Fig. 6). However, for both defects the tightest binding sites are located near the (010)  comparable, and {2H M2 } X bind as strongly to sites in the upper region as they do to sites in the lower region. dislocations, so that creation of a protonated M1 vacancy near the dislocation core is still more favorable than creation of a protonated M2 vacancy. Assuming that vacancy-related defects can lubricate glide of dislocations in olivine, it is probable that the effect will vary with the distance of the vacant site from the glide plane. In forsterite, this implies that M2 vacancies will have a greater lubrication effect for dislocations gliding on (010)  In this section, we will attempt to quantify the degree to which protonation changes segregation energies.

Comparing segregation energies for M1 and M2 defects
The degree to which the segregation energies for two defects around a particular dislocation are similar to one another can be quantified by computing a similarity measure for the segregation energy surfaces around the dislocation core. One such measure is the cosine similarity measure, which is computed for two vectors x 1 and x 2 as The similarity s 12 = -1 when the vectors are anti-correlated, while s 12 = 1 for perfectly correlated vectors. The cosine similarity measure is widely used in data mining to compare data sets, with applications ranging from facial verification (e.g. Nguyen and Bai, 2010), to comparing linguistic data sets (e.g. Liao and Xu, 2015), and automated text classification (e.g. Song et al., 2009) length equal to the number of sites, whose entries correspond to the segregation energies of each site.
Thus bare and protonated vacancies around the same dislocation can be compared provided that segregation energies have been computed for the same list of sites, as is the case in this study.
However, the similarity measure cannot be straightforwardly compared between slip systems, as the list of sites will be different. Computed values of s 12 for the M1 and M2 sites within 15 Å of the dislocation line are given in Table 3.
The similarity s 12 of the M2 sub-lattice is strictly positive for all dislocations, meaning that the  Fig. 7. Considering the segregation energies for the M1 sites at y = 0.0 (Fig. 7ab) and y = 0.5 (Fig. 7cd)

Defect segregation and olivine deformation
Vacancy-lubrication of dislocation glide has been reported in a range of different materials. Generalized stacking fault energy (GSFE) parametrized Peierls-Nabarro calculations have suggested that interstitial H may facilitate dislocation glide in Al meta (Lu et al. 2001), while the presence of interstitial O in hyper-stoichiometric UO 2 (i.e. UO 2+x , x> 0) is known to reduce the critical resolved shear stress (Keller et al. 1988), an effect attributed to interactions between the interstitial impurities and the dislocation core (Ashbee and Yust, 1982). One possible explanation is that interactions between the dislocation core and an adsorbed vacancy defect reduce the Peierls stress, although the precise mechanism remains unclear. Deformation experiments in the glide-controlled creep regime show that the critical resolved shear stress decreases from 3.8-15.0 GPa in dry olivine (Idrissi et al., 2016;Demouchy et al., 2013) to 1.6-2.9 GPa for olivine under water-saturated conditions (Katayama and Karato, 2008). This CRSS represents the stress required for deformation at 0 K, and is referred to by Katayama and Karato as the Peierls stress, although it actually represents a weighted average of the Peierls stresses for several active slip systems.
The solubility limit of vacancy-related defects in the olivine crystal lattice is relatively high, and can reach nearly 0.9 % for protonated vacancies at 12 GPa pressure (Smyth et al., 2006). However, these concentrations are probably not sufficiently great to create Peierls stress reductions of the magnitude reported by Katayama and Karato (2008). However, the strongly negative segregation energies calculated for both edge and screw dislocations mean that the concentration of vacancy-related defects will be many times greater in the dislocation core than in the bulk crystal lattice. It follows that the influence of vacancy-related defects on the deformation of olivine in dislocation-controlled creep regimes can be significant, even at low bulk concentrations. Moreover, {V M2 }″ and {2H M2 } X defects were found to have considerably lower segregation energies than equivalent defects on the M1 sublattice, and the relative abundance of M2 defects will be much higher in dislocation cores than in the bulk lattice. This is significant, as M2 vacancies are expected to have the greatest influence on the Peierls stress for dislocations gliding on (010), as these dislocations glide on the median plane of the sheet of M2O 6 octahedra.

Conclusions
Vacancy related defects are important for understanding the material properties of olivine. The addition of small quantities of water to Fo 90 olivine deforming in the glide creep regime increases strain rates, indicating a reduction of the Peierls stress. This has been plausibly attributed to lubrication of dislocation glide by protonated cation vacancies interacting with the dislocation, a process similar to the vacancy lubrication phenomenon invoked to explain flow stress variations for a range of materials.