Modelling [1 0 0] and [0 1 0] screw dislocations in MgSiO3 perovskite based on the Peierls–Nabarro–Galerkin model

In this study, we model the core structure of screw dislocations with [1 0 0] and [0 1 0] Burgers vector in MgSiO3 perovskite, in the pressure range of Earth's lower mantle (25–130 GPa). We use a generalized Peierls–Nabarro model, called Peierls–Nabarro–Galerkin, based on generalized stacking-fault energy calculations. These stacking-fault energy calculations are performed using a pairwise potential parametrization and compared to ab initio results. The results of Peierls–Nabarro–Galerkin calculations demonstrate that [1 0 0] dislocation and [0 1 0] are, respectively, characterized by a planar core spreading in (0 1 0) and (1 0 0). Our results emphasize the role of [1 0 0](0 1 0) and [0 1 0](1 0 0) slip systems in the deformation mechanism of MgSiO3 perovskite. Furthermore, we validate the use of pairwise potential for further dislocation modelling of such minerals at the atomic scale.


Introduction
Volcanism and earthquakes are spectacular manifestations of the internal activities of our planet that occur on a human timescale and that have a huge impact on mankind. Other geological Content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. phenomena like mountain building and continental drift are only noticeable on geological timescales (billions of years). Since the 1960s, the study of plate tectonics has provided a general and coherent framework to rationalize most large-scale geological observations. A more global view is necessary, however, to understand the dynamics of the Earth. Our planet is still hot: 98% of its volume is at temperatures in excess of 1000 • C and the total heat flux out of the Earth is estimated at about 44 TW. Transferring this heat through the mantle (made of thermally insulating rocks) can only be achieved by convection in the solid state. The lower mantle (between depths of about 670 km and 2900 km) is primarily constituted of a dense magnesium-rich silicate ((Mg,Fe,Al)(Si,Al)O 3 ) with a distorted perovskite structure. How this silicate creeps over geological times under mantle conditions is thus a major issue for geodynamics. In this study, we will focus on the magnesium end-member (MgSiO 3 ) of this high-pressure phase (hereafter referred to as Mg-Pv).
The pressure range of the lower mantle (25-130 GPa) represents a major challenge to study the plastic behaviour of Mg-Pv under relevant mantle conditions. Consequently, very few experimental deformation studies have been performed on this phase (Merkel et al 2003, Cordier et al 2004, Wenk et al 2004, Miyajima et al 2009. Alternatively, recent studies have demonstrated that numerical multiscale modelling can provide complementary information on plasticity of high-pressure minerals under extreme conditions . Previous numerical models of dislocations in Mg-Pv were based on a one-dimensional Peierls-Nabarro (PN) model (1D-PN: Carrez et al 2007, Ferré et al 2009. Indeed, the PN model (Peierls 1940, Nabarro 1947 represents an elegant theoretical concept for modelling dislocation core structures that has heavily influenced dislocation theory. Its attractiveness has been reinforced with the development of atomistic calculations that allow us to calculate the (non-elastic) restoring forces in the core from generalized stacking-fault (GSF) calculations, also called γ -surface (see below). The PN model has therefore been applied to a large number of crystal structures, including minerals (Joos et al 1994, von Sydow et al 1999, Lu et al 2000, Lu 2005, Miranda and Scandolo 2005, Carrez et al 2006, Durinck et al 2007. However, it remains limited by some intrinsic assumptions (Schoeck 2005). In the 1D-PN model, only planar cores (in the glide plane) can be modelled, with possible dissociation restricted to collinear configurations. Hence, the GSF is calculated along the Burgers vector only (γ -lines). The drawback of these limitations is that the model potentially predicts different solutions for screw dislocations depending on the glide plane considered, since core spreading is forced to lie in this plane (see, e.g., figure 3 and table 3 in Ferré et al 2007). The Peierls-Nabarro-Galerkin (PNG) method (Denoual 2004(Denoual , 2007 has been developed to overcome the limitations of the 1D-PN model. It is a generalization of the PN model for which non-collinear dissociations can be calculated by introducing two-dimensional GSFs (γ -surfaces) and non-planar cores can be modelled by taking into account several γsurfaces simultaneously. Using this technique, it has been shown in MgO that 1/2 1 1 0 screw dislocations in MgO change their glide planes with increasing pressure .
The goal of this study is to take advantage of the new possibilities of the PNG model to reinvestigate the structure of screw dislocations in MgSiO 3 perovskite. We focus on [1 0 0] and [0 1 0] dislocations that exhibit the lowest lattice friction on edge characters  and have already been characterized experimentally (Cordier et al 2004). We want to clarify the core structure of these screw dislocations (and hence their actual glide planes) at several pressures relevant for the Earth's lower mantle where the magnesium-rich silicate perovskite is stable. Due to the increased computational costs between γ -surfaces and γ -lines, atomistic calculations are performed here with an empirical potential (Oganov et al 2000). Since this empirical potential has never been used for calculating non-equilibrium (shear) configurations, a comparison with density functional theory (DFT) calculations is presented.

The Peierls-Nabarro-Galerkin method
In the PNG model, the dislocation core structure naturally emerges from the minimizing of an elastic energy and an interplanar potential. To illustrate the method, let us consider a unique slip plane corresponding to a potential spreading plane of a dislocation core in a volume V . Continuous deformation around the dislocation core corresponds to a three-dimensional displacement field u(r). In the same time, non-linear behaviour within the dislocation core can be represented by a two-dimensional displacement field f (r), expressed in the normal basis of the plane. Thus, the problem of core spreading consists of minimizing the energy E with respect to u and f according to equation (1), where corresponds to the material density: (1) E e corresponds to the elastic strain energy, whereas E isf is the inelastic stacking-fault energy. The latter is a function of the γ -surface energies (to be defined in the next section) from which the linear elastic part has been subtracted. Finally, minimization with respect to f is achieved by means of a time-dependent Ginzburg-Landau equation, whereas an element-free Galerkin (Zienkiewicz and Taylor 2000) method is used to compute the evolution of u(r). For the sake of clarity, the Galerkin method relies here on a two-dimensional nodal mesh (corresponding to mapping the x, y plane surrounding a dislocation aligned along the z axis). In doing so, the method allows us to also take into account several potential glide planes simultaneously and thus calculate complex (possibly three-dimensional) dislocation cores. Extended details of the PNG method are available in Denoual (2007), Pillon et al (2007) or Pillon and Denoual (2009). The nodal meshes used in this study are built with respect to the Pbnm symmetry of Mg-Pv (Horiuchi et al 1987) considering the first four low index crystallographic planes to control the periodic variation of the dislocation core energy as a function of its position (see the online supplementary materials section, supplementary figure 1 stacks.iop.org/MSMSE/22/025020/mmedia). For all calculations, mesh dimensions are equivalent to 30 perovskite unit cells with a nodal resolution of 16 nodes per Burgers vector. Once nodes supporting γ -surface energies have been attributed, a discrete dislocation is introduced into the volume. Boundary conditions consistent with a dislocation in an infinite medium are used by imposing a convolution of an elementary elastic solution with the dislocation density (see Pillon and Denoual 2009 for implementation details). The equilibrium of displacement jump field f and the density of dislocation ρ, as defined as the derivative of f by the position coordinates, are finally determined through a viscous relaxation scheme. Accordingly to Denoual (2007), we check that increasing mesh dimensions or node numbers does not influence the results. The first calculations correspond to the relaxation of the dislocation core structure. In the second step, we evaluate the Peierls stress corresponding to the relaxed core. At this stage, a homogeneous deformation is progressively applied at a velocity that allows quasi-static equilibrium so as to induce an applied shear stress on a unique given glide plane. During this loading stage, noticeable evolution of the core structure may occur prior to a strong relaxation of the measured stress. Such relaxation of the measured stress is associated with a rapid evolution of the dislocation core structure followed by a displacement of one (or more) Burgers vectors. In the following, the Peierls stress will be associated to this ultimate macroscopic stress (Metsue et al 2010).

Atomistic simulation details
Classical static simulations were performed using an interatomic pairwise potential (equation (2)) taking into account long-range and short-range interactions through Coulombic interactions and Buckingham forms, respectively. Short-range interactions include repulsive and attractive van der Waals interactions: (2) In equation (2), r ij corresponds to the distance between ions of charge z i ; b ij , ρ ij and c ij are constant parameters describing the short-range interactions. We use the parametrization proposed by Oganov et al (2000), known to be particularly robust for high-temperature and also high-pressure calculations (Oganov et al 2000, Liu et al 2007, Ito and Toriumi 2010, Chen et al 2012. This parametrization is described in table 1. Note that partial charges are used and cation-cation interactions are neglected at short range. Simulations have been performed using standard programming packages GULP (Gale 1997) and LAMMPS (Plimpton 1995). Both packages rely on Ewald summation methods (e.g. In't Veld et al 2007) for Coulombic interactions. Ground states properties of Mg-Pv (lattice parameter and elastic properties) were computed using GULP, whereas GSF energies were mostly determined with LAMMPS as discussed in the following section. First principles calculations of γ -surface energies have also been performed within the framework of the DFT using the VASP package Hafner 1993, Kresse andFurthmüller 1996). Following the earlier work of Ferré et al (2007), we used the generalized gradient approximation (GGA) to model the exchange-correlation contributions (Perdew and Wang 1992) and ultrasoft pseudopotentials (e.g. Vanderbilt 1990 or Kresse and Hafner 1994) for ionic interactions. All calculations have been performed using a single cut-off energy value of 600 eV for the plane wave expansion and a Monkhorst-Pack grid (Monkhorst and Pack 1976) scheme was used for first Brillouin zone sampling.

Generalized stacking-fault energy calculations set-up
The generalized stacking-fault energy (i.e. γ -surface) represents the energy cost per unit area incurred by the relative shear displacement f across a plane of two perfect crystal halves. Atomistic calculations of the γ -surface generally involve the use of supercells and periodic boundary conditions. After shearing the supercell, atomic relaxations are allowed along the direction perpendicular to the fault plane exclusively.
For γ -surface calculations using DFT (VASP), we used the supercells from Ferré et al (2007). Due to calculation time limitations with DFT, the number of atomic layers parallel to the shear plane is limited and the supercell must terminate on a layer with fixed atoms followed by a vacuum buffer. In that case, the size of the supercell might not be sufficient to allow full relaxation of atomic positions close to the shear plane. To overcome this limitation, we have also performed γ -surface calculations using an empirical potential (LAMMPS). In that case, we use fully periodic supercells built on lattice parameters (a 1 , a 2 , a 3 ). The supercells involve a repetition of at least eight perovskite unit cells along the normal to the slip plane. a 1 and a 2 correspond to the shortest crystallographic lattice vectors parallel to the γ -surface. When stacking faults are introduced by displacing the upper block with respect to the lower one, the lattice vector a 3 is tilted along f to keep periodicity. Energy minimization is then undertaken, keeping volume constant. This procedure has been followed for the (1 0 0), (0 1 0) and (0 0 1) γ -surfaces. The drawback is that this procedure can only be applied to γ -surfaces for which a 3 is aligned with a crystallographic axis. For a Pbnm orthorhombic structure, some γ surfaces involving shear in (1 0 1) or (0 1 1) do not meet this requirement. For (1 0 1) γ -surfaces, it is possible to use the high-index crystallographic axis [25 0 12] as an approximation of the normal direction (with a departure of 1 • from the true normal orientation). It is worth noting that it may induce some distortions compared to a perfect bulk system. For (0 1 1) γ -surfaces, keeping reasonably low distortions was not possible. We therefore used the same kind of supercells as for DFT calculations with a vacuum buffer along the direction normal to the shear plane. However, size is not a limitation and we used a neutral supercell containing at least 3280 atoms.

Results
An orthorhombic MgSiO 3 perovskite unit cell (containing 20 atoms, i.e. 4 formula unit) is first fully optimized at 30, 60, 100 and 140 GPa through static calculations performed with the GULP package. Regarding the elastic properties, the full athermal elastic constant tensors were derived from stress-strain relations (Barron andKlein 1965, Wallace 1972) implemented in GULP. As expected from the potential parametrizations used (Oganov et al 2000), the bulk properties (table 2) are found in good agreement with the available literature data, including the DFT calculations. Only C 33 and C 66 are slightly softer and C 13 and C 23 are stiffer. As a consequence, the energy coefficient for screw dislocations of Burgers vectors [1 0 0] and [0 1 0], calculated within the framework of the Stroh theory using the DisDi software (Douin et al 1986), are found to be lower in case of pairwise potential calculations and also less anisotropic compared to DFT-based evaluations.  (0 1 0). Increasing pressure to 100 GPa leads to a severe increase of γ -surface energies, with the lowest energy paths being multiplied by 1.5-2 (figure 2). Figure 3 shows the results of pairwise potential calculations for the five γ -surfaces (only those calculated at 140 GPa are presented). The agreement with the DFT calculations is shown in figure 2 along some γ -lines. Not only the shape (including asymmetries) but also the energies involved are well captured using the pairwise potential. Whatever the γ -surface (and applied pressure), the largest energy differences between the two calculations types correspond to unstable fault configurations (with an energy difference smaller than 3 J m −2 ), whereas in case of stable faults, the energy differences are reduced below 1 J m −2 . The pairwise potential parametrization is thus able to reproduce the asymmetry of energy profiles as well as their evolutions with pressure. Finally, for (1 0 1) and (0 1 1)

Dislocation core structures
Dislocation core structures are calculated with the PNG method. As an input of the PNG model, we used γ -surfaces calculated with pairwise potentials with the LAMMPS package, and elasticity tensors determined with the GULP package.  Note that both core structures exhibit a narrow core in (0 0 1) and a wider spreading in (0 1 0) or (1 0 0) (for the [1 0 0] and [0 1 0] Burgers vectors, respectively). This tendency for core spreading is indicated by a shoulder on the disregistry curve. As such a feature may indicate a tendency for dissociation, we performed a smooth fitting of the disregistry functions using a sum of arctan functions (equation (3)), whereas for the (0 0 1) plane, we rely on the canonical analytical solution (i.e. i = 0 and x i = 0 in equation (3) of Peierls 1940) f  (0 1 0) and (0 0 1) planes. Half-widths are given in Angstroms (also in reduced coordinate ζ /a in brackets where a is taken as the lattice distance in the direction of spreading). The separation distance between the two fractionals in (0 1 0)  The parameters describing the calculated relaxed cores are given in tables 3 and 4. The misfit distribution can be described as two asymmetrical fractional dislocations separated by a distance (corresponding to x 1 − x 2 in equation (3)). In case of the [1 0 0] dislocation, the separation width between the two fractionals is slightly lower than 4 Å, however, the fractionals remain strongly correlated due to the width ζ i of fractionals being greater than 1 Å. For the [0 1 0] screw dislocation, the tendency toward dissociation is more pronounced with two almost identical fractional peaks and a separation width of 4.5 Å.

Evolution of dislocation cores with pressure and Peierls stress calculations.
For calculations performed under pressures above 30 GPa, the main characteristics of the dislocation cores remain unchanged, i.e. we still find fractional decomposition of the [1 0 0] screw dislocation in (0 1 0) and of the [0 1 0] screw dislocation in (1 0 0). At the same time, core spreading in (0 0 1) remains very limited with a reduced half-width (ζ /a , where a is the lattice distance in the spreading direction) around 0.1 (tables 3 and 4). Core calculations performed at 30, 60, 100 and 140 GPa are presented in figure 6. Increasing pressure induces a noticeable decrease of the fractional width separation (at 140 GPa, fractionals of [0 1 0] dislocation are separated by less than 3 Å compared to the initial separation distance of 4.5 Å at 30 GPa). As mentioned in section 2.1, the Peierls stresses (σ p ) for the four slip systems investigated here, [1 0 0](0 1 0), [1 0 0](0 0 1), [0 1 0](1 0 0) and [0 1 0](0 0 1), have been evaluated by shearing the PNG nodal mesh. For each slip system, we applied positive or negative shear component. We observe distinct behaviours with respect to the sign of the applied stress. A typical example of core modification prior to displacement is given in figure 7. Therefore, we define here two Peierls stresses for each slip system (summarized in table 5). The difference between positive and negative Peierls stresses can be substantial, with largest discrepancies found in (0 0 1). Nevertheless, whatever the pressure or the direction of the applied strain, stresses required to move [1 0 0] and [0 1 0] dislocations in the (0 1 0) and (1 0 0) planes, respectively, are always found to be significantly lower than the stresses for glide in (0 0 1).

Discussion
Dislocation core calculations generally imply a cluster or dipole method to determine atomistic arrangement within the core (see, e.g., Woodward 2005 or Vitek andPaidar 2008). Such calculations are largely developed in material science and used for core structure calculations in metal or semiconductors (Yamaguchi and Vitek 1973, Xu and Moriarty 1998, Woodward and Rao 2002, Wang et al 2003, Pizzagalli and Beauchamp 2004, Ventelon and Willaime 2007, Groger et al 2008, Clouet 2009, Pizzagalli et al 2009, Proville et al 2012. However, for minerals and more generally ionic compounds, the number of studies is still limited (Woo and Puls 1977, Watson et al 1999, Walker et al 2004, Hirel et al 2012. A main reason is related to the technical difficulty in tracting large-scale simulations on material already containing tens of atoms in a unit cell and more than two atomic species (Walker et al 2010). Atomistic calculations of dislocation core structures in minerals or ceramics therefore rely on the use of empirical potentials. In such a context, the present study suggests that pairwise potential can be used for atomistic calculations of dislocation cores in Mg-Pv. The first point of comparison is between γ -surface calculations based on pairwise potential and DFT calculations. Whatever the pressure investigated, we find that pairwise potential calculations are able to capture the energy landscape (shape and energy level) of γ -surfaces calculated with DFT. In particular, the non-symmetric shape of γ -surfaces (figures 1 and 3) attributed to the tilting of octahedra in orthorhombic Mg-Pv (Ferré et al , 2009) is well reproduced by the pairwise potentials. The most significant discrepancies come from the maximum unstable energies along the [0 0 1] shear directions in the (1 0 0) and (0 1 0) planes (figure 2). DFT maximum energies are indeed higher without actually showing large discrepancies for the stable stacking energies at 1/2[0 0 1]. We believe that this can be attributed to a size effect of the DFT simulation cells. Indeed, the supercells used for the DFT calculations Figure 7. Illustration of the non-symmetric behaviour of [0 1 0] dislocations submitted to an applied positive or negative strain to induce glide in (1 0 0). In this example, a negative applied strain induces a first modification of the lowest fractional, which starts to move for ε − = 0.0048. It is then followed by a rapid evolution of the whole core structure with an irreversible displacement at ε − = 0.0052. On the contrary, for an applied positive strain, the largest fractional seems to be stable up to ε + = 0.0108, prior to the whole displacement of the core at 0.0112. correspond to two unit cells (40 atoms) only, including frozen atoms on supercell surfaces. This leaves less than a half unit cell to accommodate the imposed shear and relax the energy. Such supercell geometry is thus expected to overestimate unstable configuration energies. In the same time, stable stacking-fault configurations appear less sensitive to the degree of freedom allowed with a better agreement between the DFT and pairwise potential calculations. Therefore it seems that the parametrization of pairwise potentials used here is demonstrated to be able to reproduce the shear properties of the Mg-Pv structure under large strains beyond the classic elasticity limits used by Oganov et al (2000) to fit the potentials.
In particular, differentiating the Peierls stress from the misfit energy does not account for possible core modifications under stress and, accordingly, variation of the elastic energy. The PNG model accounts for both core modifications and elastic energy variations. It is also worth noticing that accurate Peierls stress calculation remains a challenging issue and may require strong refinement, as shown in a recent study by Proville et al (2012). At this stage, it is therefore preferable to consider the differences between Peierls stresses of different slip systems calculated within the same model. Here, this highlights the importance of the [1 0 0](0 1 0) and [0 1 0](1 0 0) slip systems in Mg-Pv within the pressure range of Earth mantle.
Regarding the core behaviour under stress, the stress (or strain) threshold for dislocation displacement seems to be related to the dissociation into distinct fractionals. The lowest Peierls stress value systematically corresponds to the first displacement of the widest fractional, whereas the highest value is reached when the dislocation displacement involves the first modification of the narrowest fractional (figure 7). Further investigation of this mechanism is, however, difficult to achieve by PNG as a whole displacement of the core occurs very rapidly after the first event. Therefore, further atomistic calculations are mandatory to understand this phenomenon. We believe that atomistic calculations based on the pairwise potential parametrization used here could help clarify these observations.

Concluding remarks
Based on the generalized PN model, this study clarifies the core structure of screw dislocations of Burgers vectors [1 0 0] and [0 1 0] in magnesium silicate perovskite (Mg-Pv) in a pressure range of 30-140 GPa. Screw cores are found to spread in (0 1 0) and (1 0 0), suggesting the major role of the [1 0 0](0 1 0) and [0 1 0](1 0 0) slip systems in this orthorhombic perovskite structure. In contrast to previous studies , Mainprice et al 2008, Ferré et al 2009, current calculations rely on an empirical parametrization of pairwise potentials (Oganov et al 2000). We show that the pairwise potentials calculations of γ -surface energies compare well with the DFT calculations. The accuracy of the pairwise potentials to model nonequilibrium configurations, especially the large atomistic displacements involved in dislocation cores, opens a new route to full atomistic calculations of the dislocations and mechanical properties in Mg-Pv.