The discontinuous effect of pressure on twin boundary strength in MgO

MgO makes up about 20% of the Earth’s lower mantle; hence, its rheological behaviour is important for the dynamics and evolution of the Earth. Here, we investigate the strength of twin boundaries from 0 to 120 GPa using DFT calculations together with structure prediction methods. As expected, we find that the energy barrier and critical stress for shear-coupled migration of the 310/[001] interface vary strongly with pressure. However, what is surprising is that the twin boundary also exhibits sudden strong discontinuities in strength which can both weaken and strengthen the boundary with increasing pressure. Since twin boundary migration is a proposed mechanism for both deformation and seismic attenuation in MgO, these results may suggest that MgO can undergo sudden changes in rheology due to transitions in grain boundary structure. The multiplicity of interfaces, however, necessitates the need for further studies to examine the role that phase changes in grain boundary structure play in mediating polycrystalline plasticity in the Earth.


Introduction
The Earth's lower mantle is thought to be composed of three main phases: (Mg,Fe)SiO 3 bridgmanite (Mg-pv), calcium perovskite (Ca-pv) and (Fe,Mg)O ferropericlase (Ringwood 1966). Although volumetrically inferior to the combined silicate phases, MgO is relatively ductile and so an interconnected network may control the bulk rheology of the Earth's mantle (Girard et al. 2016;Miyagi and Wenk 2016). Indeed at 20 volume per cent, MgO is close to the percolation threshold; hence, its rheology is likely to be important in areas of the mantle undergoing high strain. MgO's comparative low strength as revealed through high-pressure deformation experiments can be explained by low energy barriers required for dislocation and vacancy propagation; however, the role of grain boundaries (GB) may also be an important component, especially if the grain size is small as predicted by grain growth studies (Yamazaki et al. 1996).
This study investigates the mechanical strength of twin boundaries in MgO in relation to their potential significance as a deformation mechanism within Earth's interior. Twin boundaries are important as their mobility is a mechanism for anelasticity and seismic attenuation, but also their structure, strength and mobility are commonly used as an analogue for more general grain boundaries. Work by Harris and co-workers, as well as more recent ab initio studies on a handful of {n10} interfaces have illustrated the influence of pressure on fundamental properties such as interfacial volume V GB and energy E GB (McKenna and Shluger 2009;Verma and Karki 2010;Harris et al. 1996;Yokoi and Yoshiya 2018). These, previous studies have also shown that similar to single crystals, interfaces may, at certain P-T conditions, undergo a series of structural phase transitions.
Unconstrained by system size, a recent empirical potentials study by Hirel and co-workers has expanded the range of studied interfacial misorientations to over 200. This extensive study not only matched geometric models for lowangle grain boundaries but also, established statistically significant trends on how the effects of pressure on interfacial energies and volumes vary as a function of misorientation (Hirel et al. 2018). However, the dynamical behavioursstrength and mobility-of these interfaces have not been examined, despite their significance in the Earth's mantle.
Shear-coupled grain boundary migration, Fig. 1, has been well established as an important mechanism for both grain growth and deformation in polycrystalline assemblages (Molodov et al. 2015;Cahn and Taylor 2004;Cahn et al. 2006;Suzuki and Mishin 2005;Ivanov and Mishin 2008;Molodov and Molodov 2018). For a given shear displacement along an interface or boundary, the boundary can move normal to the shear direction allowing both grain growth and deformation to occur at the same time under stress. In some materials, this deformation mode is more effective at enhancing plastic deformation than more common deformation mechanisms such as pure grain boundary sliding (Molodov and Molodov 2018). However, it is also becoming clear that the character of the particular grain boundary or interface plays a significant role in determining the efficiency of shear-coupled migration, with different grain boundaries having very different shear strengths (Combe et al. 2017).
In this paper, we examine shear-coupled migration in MgO as a function of pressure. In order to model the dynamical processes of an interface such as a grain boundary, determining the lowest energy interfacial atomistic arrangement is required. This is non-trivial as grain boundaries are large and contain many atoms. Due to the complexity of interfaces, we use ab initio random structure searching to ensure the potential energy landscape has been thoroughly searched. Once the ground state structure of the interfaces at a given pressure is found, the Climbing-Image Nudge Elastic Band (CINEB) technique (Henkelman et al. 2000) is used to explore the migration energy barrier for shear-coupled migrations. Using these methods, we find that MgO exhibits strong discontinuities in the strength of the 310/[001] twin boundary, and over parts of the lower mantle this twin boundary is competitive with dislocation mobility, potentially allowing for fast grain growth, high attenuation and low viscosity.

Structure searching
Due to the large number of atoms required to facilitate both the periodicity of an interfacial structure and a sufficient separation between interfacial planes, the study of interfaces within the density functional framework has inherent computational cost. This computational demand is exacerbated when performing complimentary techniques such as CI-NEB calculations (Henkelman et al. 2000); thus, the scope of this study is limited to a single Coincidence Site Lattice (CSL) twin boundary, the 310/ [001]. This twin boundary has been subject to a wide range of previous studies atomistic studies as well as featuring prominently in experimental studies (Kingery 1974). Further, the 310/[001] interface has been suggested to undergo a series of structural phase transitions at discrete P-T conditions. To Fig. 1 This figure shows a simple representation of shear-coupled migration modified from Cahn and Taylor (2004). During the translational movement of length u ∥ = a n 10 [n10] , the twin boundary of misorientation 2 will migrate a distance of u ⊥ into the adjacent grain. Here, the diagonal lines represent planes of periodically repeating atoms investigate each of these successive transitions (McKenna and Shluger 2009;Verma and Karki 2010;Harris et al. 1996;Yokoi and Yoshiya 2018), we have employed the AIRSS code, which has been used to predict grain boundary structures of both graphene and strontium titanate, SrTiO 3 (Schusteritsch and Pickard 2014).
The fully ab initio method, employed by AIRSS Needs 2006, 2011), relies on placing atoms at random but 'sensible' locations within the grain boundary region that divides two MgO slabs at a given orientation. At each pressure, eight formula units of MgO were added to the free space region of height 4.0Å and cross-sectional area shown in Table 1. A minimum separation between atomic species is used to increase the efficiency of grain boundary structure searching. By imposing a minimum separation of 1.5 Å , a large number of predictably high energy structures are discarded before any ionic relaxations are performed. The remaining lowest energy structures are then relaxed in CASTEP (Hohenberg and Kohn 1964;Kohn and Sham 1965;Payne et al. 1992;Clark et al. 2005) using the Broyden-Fletcher-Goldfarb-Shanno method, such that an orthorhombic geometry is preserved. All forces are relaxed to a force tolerance of 0.05eV∕Å using an energy cutoff of 350 eV over [2 × 1 × 1] Monkhorst-Pack mesh centred on the gamma point. The PBE functional (Perdew et al. 1996;Paier et al. 2005) was used in all cases. Structural searches at a given pressure produce over 100 different GB structures from which the lowest five structures are then relaxed using VASP (Kresse and Furthmüller 1996); the code in which the CI-NEB (Henkelman et al. 2000) was performed. A full description of the above method can be found at (Schusteritsch and Pickard 2014). The nature of structure searching relies on sampling a large number of atomic arrangements. Although this imposes limitations in the precision of electronic relaxations, a 350eV cutoff was found to sufficiently distinguish a set of relaxed structures.
Electronic and ionic relaxations with VASP were performed using PBE functional (Mg pv , O) (Hohenberg and Kohn 1964;Kohn and Sham 1965;Kresse and Joubert 1996;Perdew et al. 2008), to an electronic tolerance of 10 −5 eV and force tolerance of 10 −4 eV∕Å . The wave function was solved over a 2 × 1 × 1 k-point mesh centred on the gamma point to an energy cutoff of 450 eV.
As with most atomistic studies of grain boundaries, each simulation cell used for structure searching, relaxation and CINEB calculations contains two identical parallel grain interfaces. The interfacial free energy E GB is defined as the difference in energy of a defective cell and the energy of a perfect lattice containing the same number of constitutive atoms: (1) in which E Int is the energy of the grain boundary structure, E Single is the energy of the perfect single crystal structure at equivalent pressures with equal number of constituent atoms and A Int is the cross-sectional area of the interfacial region of the defective cell. Similarly, the interfacial volumes V GB are defined as: The perfect crystalline structure was calculated with an energy cutoff and k-point sampling density, equal to that used in the defective cell.

Migration energy barrier
Shear-coupled migration is defined by a horizontal displacement of two grains which is accompanied by a perpendicular migration of the interface. This is shown schematically in Fig. 1. The periodic cells used to simulate shear-coupled migration contain two grain boundaries. During shear-coupled migration, both interfaces maintain a constant separation as they have the same sense of movement perpendicular to the shear direction. Shear-coupled migration is facilitated by the atomistic reconstruction of the grain boundary, as shown in Fig. 5. The initial and concluding states of the shear-coupled migration are determined as defined by the lowest energy structures, Table 1. An initial set of intermediate coordinates is then found by interpolating between initial and final structures, with the use of the VTST toolkit. For a given structure, the translational distance can be defined as u ∥ = a 0 10 [310] , where a 0 is the lattice vector. In which an orientation specific ideal coupling factor is defined by the translational component, u ∥ and the perpendicular migration component u ⊥ , such that: During the course of shear-coupled migration, the energy barrier ΔH is defined from the ground state energies E G and saddle point energies E S at pressure P is thus defined as Three images plus the initial and starting positions are sufficient to fully constrain the potential energy surface of the reaction. These interpolated structures act as the best initial estimate for the CINEB calculation. Here, the CINEB algorithm ensures that the migration pathway is forced over the lowest point in the potential energy landscape and by extension the lowest energy positions of the interfacial atoms. Constant volume was maintained throughout a single CINEB calculation, ensuring the PV term for each mechanism is kept constant, and thus, only differences in the enthalpy barrier of the mechanism are explored (Catlow et al. 1981). Reaction pathways were later relaxed to pressures equal to that of the initial and final images. We find only small deviations in enthalpy barrier, derived at constant pressure and at constant volume. The CINEB calculation was performed using a spring constant of −0.5 eV∕Å , in which ions were relaxed using the conjugate gradient method. During the course of the CINEB calculation, twenty atoms within the centre of the two grains are fixed, with all atoms within and around the grain boundary region being allowed to relax. As with initial structural relaxations, electronic relaxation was solved using the standard PAW-PBE functionals with an energy cutoff of 450 eV and a 2 × 1 × 1 k-point grid, centred on the gamma point. Following the convergence of the CINEB calculations, the VSPTC toolkit was further utilised to interpolate a force based cubic spline derived from extrema found along the band. This process facilitates the creation of a smoothly varying energy barrier over a set of 100 intermediate reaction coordinates. Equivalent energy barrier profiles were found when comparing to CINEB calculations conducted with five and seven intermediate images. This structure contains large voids, allowing for cationanion bonding across the interface and is stable to about 5 GPa. The interfacial energy of 1.47 J m −2 at 1.23 GPa found in this study is very similar to the previous ab initio study of 1.51 J m −2 by Verma and Karki (LDA) and somewhat below that of 1.95 J m −2 found in the embedded cluster study of McKenna and Shluger. There is also a small discrepancy between this work and the pair potential study of Hirel et al. (2018), who found an interfacial energy of 1.63 J m −2 . Despite the LDA functional  Harris et al. 1996;Yokoi and Yoshiya 2018) and previously thought to be the most stable structure at intermediate pressures. Structure (c) is more compact than the previous structure (a), with a small shift in the GB plane allowing Mg-O bonds to form whilst reducing the volume of the interface. Finally, structure (d) stable at pressure > 100 GPa is the most stable high-pressure structure and is observed in this and previous work (Verma and Karki 2010;Harris et al. 1996;Hirel et al. 2018). It features a relative shift in atomic planes parallel to the interfacial axis of rotation, resulting in a structure containing dislocation pipes. consistently producing slightly higher interfacial energy about 0.3 J m −2 , the low-pressure results are generally in good agreement with each other, with all methods predicting the same low-pressure grain boundary, and with a zero pressure interfacial energy of about 1.5 J m −2 .

Interfacial volumes and energies
At pressures greater than 5 GPa , previous studies (McKenna and Shluger 2009;Verma and Karki 2010;Harris et al. 1996;Yokoi and Yoshiya 2018;Hirel et al. 2018) suggest a transition to an atomically denser, symmetric structure, referred to in this study as structure (b) (see Fig. 2). Despite featuring within the structure search of this study, we find a lower energy structure, structure (c), to be more stable over the pressure range of 5 to 90 GPa. The differences in energy between structure (b) and structure (c) are manifested as the migration energy barrier for shear-coupled migration, as the previously theorised structure (b) defines the saddle point of the reaction. We suggest that the extensive structure search facilitated by the use of the AIRSS code allowed us to identify this lower energy interface missed in other studies. By inducing a relative shift parallel to the interface, a decrease in the volume of the structure and a continuation of Mg-O bonding across the interface can occur. Unlike previous studies, this lower energy structure avoids unfavourable cation-cation interactions across the interfacial core.
The discovery of this new structure also acts to reduce stability field of the open structure (a) to about 5 GPa , a quarter of the size suggested by previous studies. The final structural transition occurs at pressures between 90 and 100 GPa, in which a dense structure (d) featuring an offset in atomic planes parallel to the grain boundaries axis of rotation, emerges as the most stable. This offset in the atomic planes has been previously referred to as dislocation pipes and exhibits similar characteristics to that of an edge dislocation (Verma and Karki 2010). This is consistent with previous studies, as shown in Fig. 3. In summary, we find the same high and low-pressure interfacial structure as in previous work, and with similar interfacial energies (where given in previous work). However, the use of structure prediction search methods has allowed us to find a new interface stable over the pressure range of the majority of the mantle. (LDA) the local-density approximation (PP) signifies the use of pair potentials and (EC) the embedded cluster method. In general, there is good agreement between this study and previous work, with the largest disparities occuring between work carried out using pair potentials and DFT. Also included in this plot is the comparison between the structures relaxed using the LDA functional compared to the PBE functional; it is generally observed that the LDA functional derives a higher interfacial energy and a lower interfacial volume than the PBE functional, such that at pressures < 2 GPa, structures relaxed using the LDA functional offer the best correlation to previous studies.

Migration energy barrier
The effect of pressure on the migration energy barrier (MEB) (or equivalently the surface) for shear-coupled migration is shown in Fig. 4. We find that the energy barrier shows discrete and discontinuous changes associated with each structural transition. The 0-120 GPa pressure range can therefore be divided into three distinct regions, each separated with a discontinuity in the energy barrier. The atomistic motion during the migration for the 5-90 GPa pressure range is shown in Fig. 5, with the additional mechanisms for structures (a, d) found in the attached supplementary information.
The first transition from structures (a)-(c) occurs at 5 GPa and is marked by a significant reduction in MEB from 515 to 32 mJ/m 2 . This low energy barrier for shear-coupled migration in structure (c) stems from the fact that the saddle point for migration within this pressure range is actually the previously suggested structure (b). In other words, shear-coupled migration proceeds through another low energy structure and the barrier to migration is very low. At about 100 GPa, the transition to structure (d) results in an increase in MEB from 26 to 199 mJ/m 2 . The elevated MEB at pressures >100 GPa is a likely result of dislocation pipes present within the ground state structure acting as a significant barrier to interface mobility.
One notable feature of the MEB within each region is its nonlinearity with pressure. For example, structure (a) increases in energy barrier from 492 to 522 mJ/m 2 , between 2.3 and 11.4 GPa, but then decreases with pressure to 403.2 mJ/m 2 at 31.4 GPa. A similar feature is also observed in the structure (c), in which the energy barrier increases from 32 mJ/m 2 at 11.4 GPa to peak at 64 mJ/m 2 at 41.4 GPa. Then, as pressure increases further to 91.5 GPa, the energy barrier decreases to 25.7 mJ/m 2 . This nonlinearity and in some cases, a weakening with pressure, results from the relative compressibilities of the saddle point and ground state evolving at different rates with pressure. This characteristic can be inferred from the study of Yokoi et al. The three regions each of which are separated by the two discontinues that occur at ∼ 6 GPa and ∼100 GPa, respectively. These three regions correlate with the stability field of the ground states shown in Fig. 2. The stability field of the large-volumed structure, structure (a) is reduced as compared to previous studies (McKenna and Shluger 2009;Verma and Karki 2010;Harris et al. 1996;Yokoi and Yoshiya 2018). The migration energy barrier for structure (a) is defined by the highest of the energy barrier for the 310/[001] interface. Structure (a) then transforms to the denser shifted structure as shown in Fig. 2b, the lowest of the three energy barriers. The low energy barrier for migration that involves structure (c) stems for similarity in structure between the ground state and the activated state. Upon transition to structure (d), an increase in energy barrier is observed due to the presence of dislocations parallel to the interface (2018), in which the equation of state of differing grain boundary structures also varies nonlinearly with pressure.

Energies
One of the most significant results of this study is the discovery of structure (c) and the subsequent reduction in the size of the stability field of the structure (a). Such a The use of * signifies the lowest energy structure at a given pressure Pressure (GPa) Structure Lattice vectors (a, b, c)  reduction in the stability field of structure (a) appears to be consistent with denser structures observed in experimental results (Saito et al. 2013;Bean et al. 2017). Generally, the low-pressure open grain boundary structure (a) is not seen experimentally. This discrepancy between experimental and atomistic calculations has previously been explained as being due to a complex history of the experimentally grown interfaces (Hirel et al. 2018), in which denser structures formed via sliding events are shown to remain metastable. Although this study suggests that a denser ground state (structure c) may be accessible at lower pressure than previously established, it still permits a small stability range for the open structure (a). When comparing to previous atomistic studies, however, the overall trend in pressure effect on interfacial energies and volumes appears to be consistent. The small differences between the interfacial energy found in the study as compared to previous work may be due to a wide range of contributing factors: from the choice of exchange correlation functional used (McKenna and Shluger 2009;Verma and Karki 2010), the particular pair potentials (Hirel et al. 2018;Harris et al. 1996;Yokoi and Yoshiya 2018), the cell size or the method employed in searching for the ground state structure (Yokoi and Yoshiya 2018).

Migration energy barrier and critical stress
The results of this study show huge and unexpected differences in the shear-coupled migration energy for the same orientation twin grain boundary (310)/[001] as a function of pressure. These large and sudden changes are caused by different grain boundary arrangements becoming stable. These twins within MgO should feature prominently within the lower mantle, formed naturally during either growth or deformation.
As a stand alone deformation mechanism, shear-coupled migration will compete with processes such as dislocation multiplication. For comparisons to be made to these alternative processes, the critical stress for shear-coupled migration must been determined. This was calculated through the differentiation of the potential energy surface with respect to the displacement horizontal migration, i.e.

Fig. 5
The set of simulation cells defining shear-coupled migration of ground state structure (c) at 30 GPa. The CINEB sequence is formed of three intermediate images, Fig. 5 (2,3,4) along with initial and terminating structure. During shear-coupled migration, the central block of atoms, bounded by the two interfaces, is translated to the right and is indicated by the large black arrows. Atomistic rearrangement of the interfaces during shear-coupled migration is illustrated by the pattern grey atoms of the upper boundary. Whilst bottom interface shows the transition in the shape of structural units of the interface as well as the relative positions of the interface before and after a shear-coupled migration. In this regime, the saddle point of the reaction defines the previously suspected ground state interface, structure (b). The similarity in structure and energy along with the small distances covered by atoms during the cooperative shuffling events culminate in lower energy barrier regime as shown in Fig. 4 between 5 and 90 GPa. The mechanism features a shear displacement of 1.26 Å and a vertical migration of 1.26 Å , thus = 1.00 Subsequently, we find that critical stresses for shear-coupled migration which range between 0.8 and 2 GPa, within the 10-90 GPa pressure range are competitive with the 0K Peierls stress for dislocation glide in MgO . The stress and temperature dependence of shear-coupled migration can be modelled through transition state theory. Following the approach of (Ivanov and Mishin 2008), the velocity of the migrating interface can be defined as in which = 10 13 s is the vibrational frequency of the atoms, G is the average grain size of the material and k B is Boltzmann's constant. The stress-dependant migration energy barrier is defined as Due to the energy barrier's dependence on G 2 , shear-coupled migration is extremely sensitive to the grain size of the material. Such that for a 0.05 μ m grain size, we find shear-coupled migration is unable to match the velocities of dislocation glide ( 10 −21 − 10 −6 m/s) at equivalent stresses . Above this grain size, shear-coupled migration on the (310)/[001] can be classified as predominantly athermal, such the interface will only migrate at stresses close to the critical stress. Such stresses may exist within the cold interiors of subducting slabs (Garel et al. 2014;Motoki and Ballmer 2015). Therefore, pure shear-coupled migration should predominantly be seen as a deformation mechanism in materials with small grain sizes experiencing shear stresses in excess of 500 MPa. However, it is worth noting that these results are only for the (310)/[001] interface. It is definitely possible that other twin boundaries may have lower barriers to shear, thus enabling the activation of the shear-coupled migration at a wider range lower mantle conditions and grain sizes ( > 10 μm ) (Yamazaki et al. 1996). Furthermore, at large grain sizes, shear-coupled migration may be facilitated by mechanisms such as disconnection nucleation (Han et al. 2018;Rajabzadeh et al. 2013;Sun et al. 2018); however, the acute dependence of such mechanisms on the shear modulus of the material will likely be a limiting factor.
At higher pressures, a significant discontinuous increase in both shear-coupled migration energy barrier and critical stress in the transition from structure (c) to (d) is observed, i.e. a transition from a flat-shifted structure to a structure containing dislocation pipes. These dislocation pipes manifest as in-plane offsets between the atomic layers of the two grains, which means atoms within the grain boundary region have to migrate a additional distance along the axis of rotation. Furthermore, the offset in atomic planes requires atoms up to 5 Å away from the grain boundary region to undergo small rearrangements. Other interfaces such as the (210)/[001] and the (410)/[001] have been shown to exhibit dislocation pipes similar to structure (d) (Verma and Karki 2010;Harris et al. 1996). Typically occurring at pressures > 100 GPa , the formation of dislocation pipes will likely continue to be a feature of interfaces throughout the lowermost lower mantle. As shown in Fig. 4, the role which these dislocation pipes have on both the migration energy barrier and critical stress is significant. If dislocations within interfaces at lower mantle conditions are prominent, it could be argued that the mechanical role of interfaces in deformation and grain growth at lowermost lower mantle may substantially diminish.
Despite this study's sole focus on the mechanism of shearcoupled migration, the discontinuous nature of the MEB and critical stress with pressure reveals an essential consideration for other large-scale mechanisms and the broader role of which structural transitions mediate the overall rheology of a material. Equally, it may be too early to infer that discontinues in grain boundary structure are wholly responsible for observed rheological transitions in the Earth's lower mantle, notably D" (Murakami et al. 2004;Dobson et al. 2013) and 1000 km discontinuity (Rudolph et al. 2015;Jenkins et al. 2017). However, this study reveals one may expect abrupt changes in the physical properties of a bulk system stemming from the small structural transitions of interfaces.

Conclusion
As with previous studies (Verma and Karki 2010;Harris et al. 1996;Yokoi and Yoshiya 2018;Hirel et al. 2018), we find that pressure changes the ground state structure of the interface. Underpinned by the interplay between volume and bonding environment, twin boundaries opt for a relatively large volume to allow for cation-anion bonding across the interface at lower pressures. Increasing pressure then results in successive transformations into denser structures, mediated by relative displacements in adjacent atomic planes. In using the AIRSS structure prediction algorithm, implemented with full ab initio density functional theory, a means of searching the bonding environment from an unbiased perspective was performed. This confirmed previously known interfaces, as well as the discovery of a new grain boundary structure. We also show that these structural transitions cause abrupt changes in grain boundary mobility, with potential implications for bulk rheology. We suggest that if sudden transitions observed in 310/[001] features in other twin and generalised grain boundaries, the rheological behaviour of polycrystalline materials may undergo discontinuous and abrupt changes with pressure. Moreover, if multiple grain boundaries transition at equivalent pressures, such changes may act to increase or decrease the strength of the mantle with depth. There is no reason, therefore, to expect mantle rheology to vary smoothly with depth, even away from phase transitions, and sudden changes in grain boundary structure might explain some of the observed rheological transitions. We therefore suggest further studies are needed which explore the role grain boundary structure plays in defining the rheology of an assemblage, both in ferropericlase and mantle minerals more broadly.