Dislocation modelling in Mg2SiO4 forsterite: an atomic-scale study based on the THB1 potential

Knowledge of the deformation mechanisms of (Mg,Fe)2SiO4 olivine is important for the understanding of flow and seismic anisotropy in the Earth’s upper mantle. We report here a numerical modelling at the atomic scale of dislocation structures and slip system properties in Mg2SiO4 forsterite. Our study focuses on screw dislocations of [100] and [001] Burgers vectors. Computations are performed using the so-called THB1 empirical potential set for Mg2SiO4. Results of dislocation core structures highlight the primary importance of the (010) plane for [100] slip dislocations. For [001] dislocations, we confirm the occurrence of a stable narrow core that evolves into transient planar configurations to glide in (100) and (010). Such configurations suggest a locking–unlocking mechanism.


Introduction
It is nowadays well accepted that convection in the upper mantle is constrained by the plastic behaviour of olivine (Mg,Fe) 2 SiO 4 , one of its main constituents. Seismological observations, such as seismic anisotropy, can be related to the development of lattice-preferred orientations, and inform us about flow in the mantle (e.g. Panning and Romanowicz 2004). However, as pointed out by Mainprice and co-workers (e.g. Mainprice et al 2005) or more recently by Trampert and co-workers (de Wit and Trampert 2015), seismological observations remain as indirect evidence of the deformation mechanism and must be handled with care. Due to its primary importance in geology, the plastic deformation of olivine has been the subject of numerous experimental studies since the mid-seventies (Phakey et al 1972, Durham et al 1977, Evans and Goetze 1979, Gaboriaud et al 1981, Gueguen and Darot 1982, Wang et al 1988, Barber et al 2010, Mussi et al 2014. Intracrystalline plastic deformation of olivine involves two types of dislocations characterised by either [100] or [001] Burgers vectors (here and throughout the study we refer to the Pbnm description for the crystal structure). The determination of slip planes is more uncertain, as several planes have been reported in deformation experiments involving either single crystals or polycrystalline material. Under conditions of high stress and low temperature, glide of [001] dislocations is usually reported in (100) Burgers vectors usually exhibit straight screw segments (Phakey et al 1972, Durham et al 1977, Evans and Goetze 1979, Gaboriaud et al 1981, Wang et al 1988. At low temperature, [100] dislocations gliding in (001) planes are also characterised by long straight screw segments (Durham et al 1977). For temperature above 1200°C, transmission electron microscope (TEM) investigations (Durham et al 1977, Gueguen andDarot 1982) show that, in (010), [100] dislocations exhibit long straight edge segments while the screw components appear now to be much shorter and curved.
More recently, numerical modelling of dislocations in pure-magnesium end-member forsterite Mg 2 SiO 4 was carried out by Durinck et al (2007). In this work, dislocation core structures have been derived using the Peierls-Nabarro (PN) model (e.g. Peierls 1940or Joós et al 1994. Although undertaken using ab initio calculations of generalised staking fault (GSF) energy, such core structures investigation suffers from the intrinsic limitations of the PN model (Schoeck 2005). In particular, using the PN model, Durinck et al (2007) could only study planar core configurations. Using an empirical potential parameterisation for Mg 2 SiO 4 forsterite, Walker et al (2005) computed the atomic arrangement around a screw dislocation line of Burgers vector [001]. By comparison between the two approaches, it was later concluded that [001] dislocations exhibit a non-planar core (comparable to that found in bcc metals for instance), which is consistent with the occurrence of the long [001] screw segment observed in TEM images (Carrez et al 2008). Nevertheless, any conclusion regarding the dislocation slip systems in olivine based on these previous numerical modellings only is speculative.
In the following, we propose to model the core structures of dislocations with Burgers vector [100] and [001] in Mg 2 SiO 4 olivine by performing atomic-scale calculations. As further described in the following section 2, we rely on the so-called THB1 potential (Price and Parker 1987) specifically derived to model forsterite. In order to account for polarisation effects in the material, the potential is parameterised with the inclusion of a core-shell (C/S) assumption. Since molecular static calculations are performed using the large-scale atomic/ molecular massively parallel simulator (LAMMPS) (Plimpton 1995), section 2 addresses how we implemented the THB1 potential for LAMMPS simulation. The atomic core structures and Peierls stresses for the two kinds of dislocations of interest are computed according to the methodology presented in section 3.
2. Implementation of the THB1 potential within the LAMMPS package THB1 is an empirical potential developed to model Mg 2 SiO 4 forsterite and its high-pressure polymorphs (Price and Parker 1987). It includes a long-range electrostatic term described through a classical Coulombic energy term, a short range repulsive term based on a Buckingham form, a harmonic three-body term to model the directionality of Si-O bonding, and a core-shell approximation to include electronic polarizability of ions. For computational efficiency, we implemented the THB1 potential into a new inter-atomic potential library freely available for the LAMMPS package. Our implementation currently relies on the following energy decomposition: The long-range Coulombic energy is calculated by application of a Wolf summation (Wolf et al 1999), which allows calculations on periodic or non-periodic systems. The electronic polarizability is modelled in the shell model approximation by considering oxygen ions as a massive core (i.e. the nucleus) connected by a harmonic spring to a massless shell (i.e. the electron cloud). Because of this approximation, the charge of the shell remains constant and independent of the deformation state. This additional degree of freedom enables the ion to polarise in response to an electric field due to surrounding ions. The coupling between a core and its corresponding shell is described by a harmonic spring term, U Spring defined as follows: where K 1 and K 2 are the shell spring constants, and r i is the separation distance between the core and its corresponding shell. The short-ranged Buckingham potential, U Short , describes the effect of the electron cloud on nearby neighbour ions. The Buckingham potential is parameterised by three constants A ij , B ij , and C ij between shells i and j. It's functional form is defined as follows: Finally, the directionality of Si-O bonding is included in THB1 through a harmonic three-body contribution given in equation (4): where k ijk B is a derivable spring constant, θ ijk is the angle between the O-Si-O three-body term, and θ 0 is the constant tetrahedral angle. THB1 was previously derived by Price and Parker (1987), and table 1 summarises the parameterisation used in this work. Within this implementation, the lattice parameters of forsterite are a=4.7874 Å, b=10.2717 Å, and c=6.0227 Å (as shown in figure 1). We also computed the nine elastic constants of forsterite, following the description of Beckstein et al (2001) for orthorhombic crystal structure. The calculated C ij are listed in table 2. It is worth pointing out that, compared to previously reported elastic constants, the Wolf summation method leads to identical results. Independent of the long-range summation techniques Buckingham potential Harmonic three-body term Table 2. Elastic constant C ij tensor (in GPa) using either Wolf summation method (this study) or Ewald summation (Gale 1997). First-principles calculations (Durinck et al 2005) are given for comparison with the THB1 potential.  (Ewald or Wolf), computed elastic constants are in perfect agreement with differences below 1%. Nevertheless, as commonly reported, empirical potential calculations do not fully render the first-principles elastic constants set, but correspond to the best compromise for forsterite (Price and Parker 1987). To further model dislocation structures at the atomic scale, one can also investigate the response of the potential to rather large and non-elastic displacements. This is commonly done by computing GSF. Therefore, prior to dislocation modelling, we perform GSF calculations and compare THB1 γ-energies with those computed from firstprinciples by Durinck et al (2005). For the GSF calculations, we follow the full periodic boundary conditions method described by Gouriet et al (2014) with Si atoms allowed to move perpendicularly to the stacking fault plane and, in order to increase the degrees of freedom of the system, with O and Mg atoms freely allowed to relax in any direction. Figure 2 shows the (010) γ-surface energy. Low-energy paths correspond to [100] and [001] (figures 2(b) and (c)), as expected from density function theory (DFT) calculations. The shapes of the GSF are also consistent between the two methods, and one may notice that γ-surface energies computed with the THB1 potential are systematically slightly higher. The analytical expressions of force fields implemented to the code and algorithm are given in the appendix.

Dislocation core modelling methods
The atomic systems used are designed to study dislocation core properties, and dislocation motion in different glide planes. We thus rely on a periodic system, defined as follows: the dislocation line is along z, and the normal to the glide plane is along y (see figure 3). Consequently, the motion of screw dislocation is along the x direction of the atomic system. Dealing with straight screw dislocation lines, the system is naturally periodic along the z direction (with a length equal to the magnitude of the Burgers vector, b, of the simulated dislocation). The periodicity along x is maintained by adding a tilt component of magnitude b/2 to this axis (Hirel et al 2014). For each investigated slip system, the Mg 2 SiO 4 crystal structure is rotated to fulfil these two conditions. However, dealing with an orthorhombic symmetry material, the y direction (i.e. the normal to the glide plane) is not necessarily aligned with a crystallographic direction [uvw]. Therefore, the atoms located at the top and bottom regions of the supercell along the y axis are kept fixed to mimic an infinite perfect crystal. With such a setup, the systems investigated can be considered as periodic slabs. To further minimise spurious elastic interaction between periodic replicas, the system sizes are gradually increased to typically 480 Å between dislocation replicas, and 120 Å between the glide plane and artificially frozen regions of atoms. This corresponds to system sizes between 13 440 atoms and 53 760 atoms. Straight screw dislocations are finally introduced into the system thanks to the classical displacement field either relying on isotropic elasticity solutions or anisotropic elasticity as implemented in Atomsk (Hirel 2015). All as-built systems are then optimised using a conjugate gradient minimisation algorithm with a stopping tolerance force set at 10 −12 eV/Å. Following the minimisation procedure, the dislocation core structures are analysed using differential displacement maps (Vítek 1968) or through the computation of atomic disregistries across glide planes.
To determine the dislocations mobility, a simple shearing is applied to the system, which results in a resolved shear stress in the glide plane. For screw dislocations this can be achieved by applying shear strain increments of 1%, ensuring quasi-static loading. At each increment, the system is relaxed using a conjugate gradient method. Consequently, the aforementioned loading conditions allow for determining a critical strain, and thus a critical stress above which the system does not respond purely elastically. The dislocation starts moving when an inflection on the stress-strain behaviour is obtained. In the following, various cell sizes (along the x and y directions) have been considered to ensure system size effects. Using molecular statics, the computed critical shear stresses are the Peierls stresses for the considered slip system. By repeating such a process for various planes, the methodology can be applied to  Figure 4 shows the differential displacements map around a [100] screw dislocation core. Based on neighbours' cations analysis, the core of the [100] screw dislocation shows a clear tendency of spreading in (010) with a core centre located between SiO 4 tetrahedra at the M2 site (we used the classical distinction M1, M2 of cation sites in olivine group structure; see for instance Deer et al 1982). To investigate the spreading of the core, one can extract the disregistry functions within (010) from the atomic positions. We observe (figure 4) a dissociation of the core into two strongly correlated Burgers vector density peaks. The distance between these peaks, ∼7 Å (table 3), is close to one lattice parameter in the [001] direction. Details of the atomic arrangements around the line indicate that the core of [100] dislocations involves some local edge displacements in the vicinity of (010) (figure 5). These local displacements do not produce any edge component of the dislocation Burgers vector but can be associated with a dilatation state of the core. Following Sun et al (2016), we use an atomistic-to-continuum cross over method to analyse the details of local edge displacements. The unrelaxed and relaxed configurations obtained during the atomistic calculations are used as the reference state and deformed state, respectively. The components of the displacement vector field, u, in the vicinity of the dislocation core are calculated for the three different sub-  (010), and correspond to outward displacements of atoms. Thus, from the computed gradient of the displacement vector, u, one can compute the corresponding strain tensor giving access to the local dilatation taken as the trace of the strain tensor. As shown in figure 6, the core of the [100] screw dislocation is associated with positive values of the trace of the strain tensor, confirming the dilatation state of the dislocation core. As mentioned in the methods section, a simple shearing was applied to the stable core to determine at which critical stress the dislocation starts moving in a given glide plane. We find two different critical stresses for glide in (010) and (001). For [100](010), the computed Peierls stress is 3.1 GPa. By contrast, when the [100] dislocation is forced to glide into (001) we find a higher value of the Peierls stress around 7 GPa. In the case of slip in (010), the [100] core structure evolves with stress. While it still spreads in (010), the core becomes narrower and does not show any evidence of dissociation ( figure 7 and table 3). Computing the energy difference between the two [100] screw configurations, it turns out that the undissociated dislocation core corresponds to a metastable state with a higher energy of 0.4 eV b −1 with Table 3. Dislocation disregistry analysis. Half widths ζ of the dislocation core are given according to a fit of the atomic disregistries with the following expressions: S(x)=b/2 +b/π.arctan(x/ζ) for a compact core, and S(x)=b/2+b 1 /π.arctan(x−x 1 /ζ 1 )− b 2 / π.arctan(x−x 2 /ζ 2 ) for a dissociated core. In the case of a dissociated core, b i , x i , and ζ i correspond to the partial Burgers vectors, the position of the partial cores, and their widths, respectively.

Screw dislocations with [001] Burgers vector
Inserting a [001] dislocation into the structure leads to a stable core configuration located between four SiO 4 tetrahedra with a dislocation line centred on a Mg cation site ( figure 8). As previously shown, the structure of the [001] stable core can be analysed with the differential displacement map. All significant displacement arrows around the dislocation line are restricted to a small area within one unit cell ( figure 8). Thus, the core appeared to be compact with no particular evidence of spreading in any preferential plane. This is confirmed by computing the disregistry functions in several planes (100), (010), or {110} (figure 9). Independent of the plane, the core profile is narrow with a width smaller than 5 Å (table 3), preserving the shape of the SiO 4 tetrahedra. We now investigate the critical stress needed to promote glide in these planes. We found that above a critical stress of 7.2 GPa, the dislocation moves over a large distance (with respect to the simulation cell size) in (010). Above 7.2 GPa of resolved shear stress, an analysis of the conjugate gradient minimiser steps shows that once the dislocation is set into motion, it exhibits a different core structure (figure 10), entirely planar, spread in (010). It is worth noting that such a planar configuration cannot be stabilised without applying stress in our simulations. Indeed, we verified that without stress a planar configuration always recombines into the compact stable core configuration (figure 8). However, once this transient state ( figure 10) is reached, the dislocation core can travel over several lattice repeats before finally falling again into another stable core configuration. For a slip in (100), we found an onset for core motion corresponding to a critical stress of 6.9 GPa. Once again, the transient states of the core (figure 11) computed according to the conjugate gradient algorithm show that the core tends to adopt a planar configuration in (100) during motion. As for the other case, the transient planar configuration in (100) cannot be stabilised without stress. Finally, for the [001](110) slip, whereas we apply a simple shear component to maximise the resolved   shear stress in (110), we find that above 7.8 GPa, the dislocation moves in (100) with a behaviour similar to that previously described. Indeed, due to the close orientation between (110) and (100), a simple shearing of the simulation cell cannot avoid inducing a significant resolved shear stress in (100).

Discussion
Our calculations of the atomic arrangements around the screw dislocation core of dislocations with the [001] Burgers vector show a compact core structure with no specific or preferential spreading in any given plane. As expected, differential displacement maps around the [001] core show typical features consistent with previous attempts to model screw dislocation cores Nevertheless, the choice made for system set-ups are different. In the study by Walker et al (2005), computations were performed on a small system size with, consequently, a potentially strong influence on rigid periodic boundary conditions. In the present work, the investigated systems are much larger and the boundary conditions applied are supposed to be less restrictive (Woodward 2005, Walker et al 2010. Also, as discussed in the methods section, size effects were taken into account by gradually increasing the system size. One of the main achievements of the present study is to be able to compute Peierls stresses for glide in (100),  stresses of, respectively, 6.9 GPa and 7.2 GPa. For [001](110), the determination of the Peierls stress is more ambiguous, but our results suggest that in this plane, the Peierls stress should be above 7.8 GPa. Whatever the glide plane, lattice friction opposed to [001] dislocation is therefore very high. The occurrence of high lattice friction can be viewed as a consequence of the compact core structure of [001] dislocation, and is therefore consistent with the experimental observations of long straight screw dislocation lines in deformed olivine samples (e.g. Kohlstedt et al 1976or Gaboriaud et al 1981. One may also notice that whatever the glide plane, our results suggest that glide in (100), (010), or (110) may be activated at comparable stress levels. Indeed, glide of [001] dislocations in (100) and also possibly in {110} have been reported from TEM observations (see for instance Gaboriaud et al 1981) or using decoration techniques (Kohlstedt et al 1976). More recently, based on electron tomography techniques, the three possible slip planes, i.e. (100), (010), and {110}, have been reported simultaneously (Mussi et al 2014), but with a greater occurrence of {110}, which is not explained by our calculations of lattice friction. At stresses higher than the Peierls stress, our calculations show that the [001] core adopts some planar transient state (figures 10 and 11). Since it is impossible to preserve the atomic configurations of these transient planar configurations when minimising the energy of the system if no external strain is applied, it is suggested that the planar core configurations for [001] dislocation are highly unstable. Since glide may involve high-energy configurations, we confirm a clear tendency of the [001] core to glide through a so-called locking-unlocking mechanism (Couret and Caillard 1989), as recently described at the atomic scale in hpc metals (Clouet et al 2015). This suggestion was initially proposed by Carrez et al (2008) based on unrelaxed core energy calculations and the PN model of core spreading. As for glide in (100) and (010), one can speculate that glide in (110) would involve a planar configuration. Although a comparison with PN model calculations may not be straightforward, it is however interesting to note that the transient core of the [001] screw dislocation in (010) does not show any evidence of dissociation whereas that in (100) indicates a tendency to dissociate (see figures 10(b) and 11(b)). Qualitatively, these two transient cores are similar to those computed using the PN model (Durinck et al 2007). This is another argument to suggest that glide in (110) would involve a planar dissociated transient core, as predicted by the PN model.
For [100] screw dislocations, our results show that the core dissociates in (010) into two collinear partial dislocations. This tendency to dissociation has never been reported; however, based on hcp lattice considerations, possible dissociations of dislocations in olivine were addressed by Poirier (1975). However, the computed width of dissociation is rather small with a distance between partial dislocations slightly larger than a lattice repeat along the [001] direction. Such a dissociation scheme can be viewed as a direct consequence of the occurrence of a stable stacking fault energy at ½[100] in the GSF calculations ( figure 2(b)). Unfortunately, such a feature was not clearly observed on γ-energy computed from firstprinciples by Durinck et al (2005), and may be a function of the relaxation procedure used during the GSF calculation. Indeed, for DFT calculations, Durinck et al (2005) computed the GSF with non-periodic atomic systems, for which surfaces induce discontinuities and potentially oscillations in the γ-surface energies depending on cell sizes. As a consequence, the dissociation configurations of the [100] core may not be captured by the PN model used by Durinck et al (2007). More importantly, the THB1 potential calculations of the [100] dislocation core show local edge displacement of atoms within the vicinity of the core perpendicular to (010), which is in agreement with displacements of atoms outward from (010) observed by Durinck et al (2005) along the γ-energy path. Therefore, the present calculations confirm the occurrence of the dilatation state within the [100] screw core. In the context of the upper Earth's mantle (where olivine is one of the most important minerals until a depth of ca. 410 km, where the pressure reaches 13 GPa), the dilatation state of the core may have strong implications on the activation of the [100] slip with respect to an increase in pressure. Indeed, as proposed by Durinck et al (2005), one effect of pressure is to work against this dilatation state, leading to an intrinsic hardening of the glide properties of [100](010).

Conclusion
Molecular static calculations relying on the THB1 potential have been undertaken to study dislocation core structures and slip in Mg 2 SiO 4 forsterite. We analyse the dislocation core structures of screw dislocations with [100] and [001] Burgers vectors. For [100] dislocations, our results highlight the importance of the [100](010) slip system resulting from a clear spreading of the dislocation core in this plane. By contrast, the [001] screw dislocation is characterised by a compact core. Moreover, glide analysis of [001] reveals the occurrence of several planar core configurations depending on the slip planes. It should be emphasised that such configurations can be viewed as higher energy transient states of the core related to the dislocation motion typical of a locking-unlocking mechanism. The calculation of the energy of the transient planar configurations is out of the scope of this study, but could be addressed in the future. As a final remark, we believe that the methodology presented in this study can be used to further investigate the effect of pressure on the core structure in olivine, and further revisit the unsolved question of the change in main slip system activities from [100] to [001] in olivine, as reported over the last decades.