# Ferroelasticity in MgSiO3 bridgmanite: Assessing the intrinsic structure and motion of (-110) twin walls with atomic-scale simulations

Pierre Hirel, Philippe Carrez, Patrick Cordier

► To cite this version:

Pierre Hirel, Philippe Carrez, Patrick Cordier. Ferroelasticity in MgSiO3 bridgmanite: Assessing the intrinsic structure and motion of (-110) twin walls with atomic-scale simulations. Scripta Materialia, 2020, Scripta Materialia, 188, pp.102-106. 10.1016/j.scriptamat.2020.07.016. hal-02921984

HAL Id: hal-02921984

https://hal.univ-lille.fr/hal-02921984

Submitted on 25 Aug 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Ferroelasticity in MgSiO₃ bridgmanite: Assessing the intrinsic structure and motion of (110) twin walls with atomic-scale simulations

Pierre Hirel<sup>a,*</sup>, Philippe Carrez<sup>a</sup>, Patrick Cordier<sup>a,b</sup>

<figure><img src="image_3.png" /></figure>

## ARTICLE INFO

Article history:

Received 16 April 2020

Revised 19 June 2020

Accepted 5 July 2020

Keywords:

Simulation

Mechanical properties

Crystal defects

Silicates

Twin boundaries

## ABSTRACT

Bridgmanite MgSiO₃ with perovskite structure exhibits a ferroelastic behaviour via the motion of {110} twin walls. We use atomic-scale simulations to investigate the properties of (110) twin walls in bridgmanite. We resolve their atomic structure and obtain the minimum energy path for their motion in a wide pressure range relevant to the Earth's lower mantle.

© 2020 Acta Materialia Inc. Published by Elsevier Ltd.

This is an open access article under the CC BY license. (http://creativecommons.org/licenses/by/4.0/)

Ferroelasticity, i.e. the ability for a material to exhibit a hysteresis in its stress-strain curve under cyclic loading, is commonly observed in metals with martensitic transformations [1–4]. It also occurs in ceramics with distorted perovskite-type structures, e.g. BaTiO₃ [5], CaTiO₃ [6–9], BiFeO₃ [10,11], LaAlO₃ [12,13], or LaCoO₃ [14,15]. The latter are ferroelectric or have octahedral distortions that break the symmetry between the different (100) directions, making possible the occurrence of twin walls (TW). Those TWs move reversibly under cyclic loading, causing ferroelastic behaviour.

Bridgmanite (Mg,Fe,Al)(Si,Al)O₃ with the perovskite structure is by far the dominant phase of the Earth's lower mantle, representing up to 80% of its composition [16]. As such, it is bound to play a major role in mantle convection and in seismic wave propagation. While the bulk properties are well constrained [17–19], the role of defects on the mechanical response remains largely unexplored and very challenging to characterize. A distinctive feature of bridgmanite is that it is ferroelastic due to its octahedral distortion [20,21]. Direct experiments are notoriously difficult because this phase is stable only under very high-pressure conditions typical of the lower mantle (26–125 GPa). Numerical modelling is also rendered difficult by the fact that it is a complex ternary phase with mixed ionic-covalent bonding, strong octahedral distortions, and other properties that depend on pressure. Because of these major issues, for a long time researchers fell back to analogue perovskite

materials that can be studied in more details in ambient conditions, like SrTiO₃ [22] or BaTiO₃ [23], hoping to extend their conclusions to Earth mantle perovskites. However, the mobility of TWs being strongly dependent on the chemical composition and bonding, it is difficult to obtain relevant quantitative information about bridgmanite by extrapolating data obtained in analogue materials.

In this work, we use atomic-scale simulations to study (110) twin walls (TW) in MgSiO₃ perovskite, the magnesium-rich end-member of bridgmanite. We characterize their atomic structure, formation energy, and mobility mechanisms, in the pressure range 30–120 GPa relevant to the Earth's lower mantle. We model MgSiO₃ perovskite in its *Pbnm* representation, i.e. the shortest lattice vector is **a** = [100], followed by **b** = [010], and **c** = [001] is the third shortest lattice vector. Lattice parameters depend on pressure, and range from about *a* = 4.65 Å, *b* = 4.76 Å, *c* = 6.72 Å at 30 GPa, to *a* = 4.38 Å, *b* = 4.53 Å, *c* = 6.34 Å at 120 GPa [19,24]. Interactions between Mg, Si and oxygen ions are described with a rigid-ion potential parameterized by Alfredsson et al. which was fitted to DFT calculations [25]. A comparison of 27 potentials by Chen et al. [26] demonstrated that this potential is one of those that reproduce the most accurately the lattice parameters and elastic moduli of MgSiO₂ in the pressure range investigated hereafter. In addition, our own calculations of generalized stacking faults show that planar defects energies are also on par with DFT calculations within an accuracy of 5%, so this potential is a good candidate for modelling large systems containing twin walls.

The lattice vectors [100] and [010] differ in length hence the directions [110] and [110] are not orthogonal, which is the reason

<footnote><sup>a</sup> Univ. Lille, CNRS, INRAE, Centrale Lille, UMR 8207 – UMET – Unité Matériaux et Transformations, Lille F-59000, France</footnote>

<footnote><sup>b</sup> Institut Universitaire de France, 1 rue Descartes, Paris F-75005, France</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: pierre.hirel@univ-lille.fr (P. Hirel).</footnote>
<figure><img src="image_4.png" /><figcaption>**Fig. 1.** Atomic structure of the (110) twin wall in bridgmanite at 30 GPa. Crystal directions on each side of the TW are indicated: in the top crystal, the directions **a** = [100] and **b** = [010] are swapped with respect to the bottom crystal. The direction **c** = [001] is normal to the figure. The left graph shows the tilt of SiO₆ octahedra about the [001] direction, and the middle graph the distance between consecutive Si atoms along the normal to the TW. The extension of the TW is delimited by dashed lines.</figcaption></figure>

for the occurrence of TWs. A TW is constructed by stacking two crystals of MgSiO₃ where the directions **a** = [100] and **b** = [010] are exchanged, as illustrated in Fig. 1. Using 3-D periodic boundary conditions, equivalent TWs form at the extremities of the simulation box. The cell size (530 Å length, 3200 atoms) was verified to give good convergence of the results. The energy is minimized with the conjugate-gradients algorithm. Isostatic pressure is applied by rescaling appropriately the lattice constants according to the equation of state of bridgmanite as described by the potential. Atomic systems are constructed with Atomsk [27], simulations are performed with LAMMPS [28], and visualization with VESTA [29].

Fig. 1 shows the system after relaxation at the imposed pressure of 30 GPa. Atomic planes of the two crystals meet with an angle of 182.66°. The angle of tilt of the SiO₆ octahedra changes sign when crossing the wall, as reported in the graph in the same figure. The width of the TW, obtained by fitting an arctan function, is 2w = 19.88 Å. The interplanar distance, measured as the distance between consecutive Si atoms normally to the TW, is increased at the edges of the TW by about 1% compared to the bulk. Note that there is no discontinuity in the local environment of atoms: silicon ions always have six oxygen neighbours, and Mg ions have twelve. The TW is characterized mainly by an inversion of the tilt of the SiO₆ octahedra. This structure is very similar to the one observed by Van Aert et al. with high-resolution transmission electron microscopy in CaTiO₃ [9].

The formation energy per surface area is computed as:

$$ \gamma_{\text{TW}} = \frac{E_1 - E_0}{2S} \qquad (1) $$

where E₁ is the total energy of the system containing TWs, E₀ that of a defect-free supercell containing the same number of atoms, and S the surface area of the TW. The factor 2 accounts for the fact that there are two TWs in the system. We obtain γTW = 1.37 J/m² at 30 GPa. Although it is quite high compared to TWs in other perovskites (e.g. 0.32 J/m² in CaTiO₃ [6]), one must bear in mind that it is evaluated under high pressure, and it remains lower than that of other planar defects in MgSiO₃, which typically range from about 2 to 7 J/m² at 30 GPa [30,31].

In order to test the mechanical behaviour of the sample containing TW, we submit it to cyclic shear strain εxy by increments of 0.01%, followed by relaxation. The system's internal stress is monitored, and the resulting stress-strain curve is reported in Fig. 2.

During the initial loading (εxy > 0), a TW first moves under an applied strain of about 0.04%, as shown in the inset of Fig. 2, corresponding to a stress about 120 MPa. This motion occurs by the rotation of octahedra, so that the plane of inversion is displaced normally to the TW by one interplanar distance (about 3.32 Å). The actual ions displacements remain very small, the largest displacement (an oxygen ion) being about 0.41 Å. Magnesium and silicon ions also relax inside their respective sites, and their combined motions result in an inversion of [100] and [010] lattice vectors.

As the strain continues to increase, the stress in the upper crystal becomes increasingly larger than in the lower one. In order to release stress, the TW at the center of the cell moves upwards and the one at the edges of the cell moves downwards, resulting in the growth of the lower crystal and the shrinkage of the upper one. The TW motion is analogous to that of ferromagnetic domain walls, caused by the rotation of magnetic moments, favouring the growth of some domains and the shrinkage of others.

Due to the finite size of the simulation cell, increasing strain eventually results in the TWs meeting each other and annihilating, leaving only a single domain of MgSiO₃. To avoid that, the applied strain is reversed at εxy = 4%, before the TWs meet, and then decreased, causing the two TWs to move in directions opposite to the previous ones and the stress to decrease. When the system reaches a stress-free state, the TWs do not return to their initial position, resulting in a spontaneous strain of about 1.85% (point S₁ in Fig. 2). Strain is decreased further down to -4%, and increased again to 4%. This cyclic loading is repeated twice, leading to the hysteresis curve reported in Fig. 2. Note that the hysteresis obtained here cannot be expected to reflect the actual behaviour of an experimental MgSiO₃ sample. Indeed, our numerical simulation contains only two TWs very close to one another and does not include temperature effects, while an experimental sample can typically contain several grains and various defects, and is deformed at finite temperature, thus producing very different stress-strain
<figure><img src="image_5.png" /><figcaption>Fig. 2. Simulated stress–strain curve obtained during cycling loading at imposed strain of MgSiO₃ perovskite supercell containing two twin walls, at 30 GPa and without accounting for temperature effects. Initially, the system is free of internal stress (point A). During loading, the TWs move to new metastable positions, so that removing the strain results in non-zero internal stress (point B). The operation can be repeated, showing the reversibility of the motion of TWs. S₁ and S₂ mark the spontaneous strain remaining in stress-free configurations. Inset: zoom in on the initial loading, showing the stress release associated with the first displacement of a TW at about 120 MPa (arrow).</figcaption></figure>

curves. Nonetheless, under certain conditions a MgSiO₃ sample is expected to demonstrate a hysteresis during cyclic deformation.

The previous imposed-strain simulation provides a first estimate of the critical shear stress for TW motion of about 120 MPa. However this estimation is rather crude, because the finite strain increment (0.01%) results in a large uncertainty on the measured stress (± 30 MPa). In order to obtain a more accurate value, we focus in more detail on the elementary mechanism, i.e. when a TW moves from one equilibrium position to the next. It is associated with an activation energy, which vanishes under a critical stress to be determined.

To determine the minimum energy path (MEP), we rely on the nudged elastic band (NEB) method [32,33]. Two configurations are extracted from the previous deformation run, corresponding to two consecutive stable positions of the TW, and are brought to a strain-free state and relaxed. Then, linear interpolation is used to construct seven intermediate images between the initial and final states. Finally, the NEB method is used to relax these images along the MEP. Fig. 3 shows the final MEP obtained at the pressure of 30 GPa, and the activation energy is about $E_a = 19.44$ meV or 7 mJ/m².

The critical stress for TW motion is obtained by computing the maximum derivative of the energy curve with respect to the TW position. At $P = 30$ GPa, we obtain a critical stress about $\sigma_c = 64$ MPa. This value is lower than the one obtained from the stress–strain curve (120 MPa, see Fig. 2), but it is expected to be more accurate because it is not evaluated under applied load, and was computed from a refined energy path.

We reassessed these properties under applied isostatic pressures of 60, 90, and 120 GPa. While the tilt of SiO₆ octahedra increases with pressure, it has little effect on the atomic structure and width of TWs. The formation energy however, does change significantly: we obtain $\gamma_{TW} = 1.94$ J/m² at 60 GPa, 2.48 J/m² at 90 GPa, and 2.94 J/m² at 120 GPa, i.e. it increases by a factor of 2 across the pressure range investigated. This seems to indicate that TW formation becomes increasingly difficult at larger depths.

Pressure also has a significant effect on the mobility of TWs. Using the NEB method as before, we find $E_a = 29.16$ meV at 60 GPa, 45.63 meV at 90 GPa, and 65.50 meV at 120 GPa. This represents an increase by a factor of 3.4 across the mantle. This evolution is reported as circles on the left-hand side of Fig. 4. We fitted our simulation data with a quadratic function of the form:

$$E_a(\text{eV}) = aP^2 + bP + c \quad (2)$$

The fitting procedure yields the following values of the fitting variables: $a = 2.82 \ 10^{-24}$ eV.Pa⁻², $b = 9.26 \ 10^{-14}$ eV.Pa⁻¹, and $c = 1.4 \ 10^{-2}$ eV. The corresponding curve is plotted on the left-hand side of Fig. 4. Because the domain of stability of MgSiO₃ is 26 – 125 GPa in the Earth's mantle conditions, extrapolation out of this range is not relevant.

Similarly, the critical stress for TW motion increases with pressure: $\sigma_c = 121$ MPa at 60 GPa, 260 MPa at 90 GPa, and 424 MPa at 120 GPa, i.e. it increases by a factor of 6.6 across the mantle. Again we fitted these values with a quadratic function:

$$\sigma_c(\text{Pa}) = dP^2 + eP + f \quad (3)$$

After fitting, the constant values are $d = 2.97 \ 10^{-14}$ Pa⁻¹, $e = -3.95 \ 10^{-4}$, and $f = 4.625 \ 10^7$ Pa. The data from the simulation and the curve associated with the fitted function are plotted as open squares in Fig. 4. These parameters for TW motion determined from our atomic-scale simulations can be useful as input parameters for coarse-grained models, like the ones based on the Landau theory that have been widely applied to TW dynamics in the past [8,35].

The presence of TWs in bridgmanite in the Earth's lower mantle may add an important contribution to the attenuation of seismic waves, as pointed out by Harrison and Redfern [12]. Indeed, TWs are expected to move at the passage of seismic waves, absorbing part of their energy in the process. Assuming that their motion follows an Arrhenius law with an activation energy $E_a$, it is meaningful to consider the dimensionless ratio $E_a/kT$, where $k$ is Boltzmann's constant and $T$ the temperature. We estimate it using
<figure><img src="image_6.png" /><figcaption>**Fig. 3.** Minimum energy path for the motion of a TW from a stable position to the next in MgSiO₃ perovskite at 30 GPa, as determined with the NEB method (see text). Empty circles are the results from the calculation; the continuous line is an interpolation with a cubic function.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Fig. 4.** Evolution with pressure of the activation energy (red circles) and critical stress (blue squares) for the motion of (110) twin wall in MgSiO₃ perovskite. The open symbols are the data obtained from atomistic simulations, and the continuous curves are quadratic functions fitted to the simulation data. Fitted functions are plotted only in the pressure range relevant to the Earth's lower mantle (26 - 130 GPa). The right panel shows the evolution of the dimensionless ratio *E*<sub>a</sub>/kT as function of depth (see text). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

the profile of pressure and temperature as function of depth from the PREM model [34]. At all depths we find that it ranges from 0.1 to 0.3 (right panel in Fig. 4), meaning that the activation energy is much smaller than thermal energy, and that TW motion is very easy and independent of temperature in the mantle conditions.

A limitation of our study is that it does not account for the effects of temperature. At finite temperature, TW motion may depend on the nucleation of kinks or terraces, with activation energies and temperature dependence that may be different from those presented here. Another important aspect is the interaction of TWs with other defects. In particular, vacancies are known to have a strong effect on the pinning of TWs [8]. Since bridgmanite is submitted to high temperatures in the Earth's mantle, vacancies are expected to play a major role in the mobility of TWs. Ferroelastic materials are also reputed to form complex intersecting domain patterns [5,35,36]. The interactions between intersecting TWs would have to be explicitely modelled to obtain a more realistic description of their role on seismic attenuation and the deformation of bridgmanite.

**Declaration of Competing Interest**

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

**Acknowledgements**

This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme under grant agreement No 787198 - TimeMan. Computational resources were provided by the DSI at University of Lille.

**References**

[1] V.K. Wadhawan, Bull. Mater. Sci. 6 (1984) 733–753, doi:10.1007/BF02744001.

[2] M. Ahlers, J. Pelegrina, Mater. Sci. Eng. A 356 (2003) 298–315, doi:10.1016/S0921-5093(03)00144-8.

[3] S. Sarkar, X. Ren, K. Otsuka, Phys. Rev. Lett. 95 (2005) 205702, doi:10.1103/PhysRevLett.95.205702.
[4] Y. Liu, D. Favier, L. Orgeas, J. Intell. Mater. Syst. Struct. 17 (2006) 1121–1126, doi:10.1177/1045389X06065237.

[5] G. Arlt, J. Mater. Sci. 25 (1990) 2655–2666, doi:10.1007/BF00584864.

[6] M. Calleja, M.T. Dove, E.K.H. Salje, J. Phys. Condens. Matter 15 (2003) 2301–2307, doi:10.1088/0953-8984/15/14/305.

[7] R.J. Harrison, S.A. Redfern, J. Street, Am. Mineral. 88 (2003) 574–582, doi:10.2138/am-2003-0411.

[8] W.T. Lee, E.K.H. Salje, L. Goncalves-Ferreira, M. Daraktchiev, U. Bismayer, Phys. Rev. B – Condens. Matter Mater. Phys. 73 (2006) 3–7, doi:10.1103/PhysRevB.73.214110.

[9] S. Van Aert, S. Turner, R. Delville, D. Schryvers, G. Van Tendeloo, E.K. Salje, Adv. Mater. 24 (2012) 523–527, doi:10.1002/adma.201103717.

[10] F. Kubel, H. Schmid, Acta Crystallogr. Sect. B Struct. Sci. 46 (1990) 698–702, doi:10.1107/S0108768190006887.

[11] S.H. Baek, H.W. Jang, C.M. Folkman, Y.L. Li, B. Winchester, J.X. Zhang, Q. He, Y.H. Chu, C.T. Nelson, M.S. Rzchowski, X.Q. Pan, R. Ramesh, L.Q. Chen, C.B. Eom, Nat. Mater. 9 (2010) 309–314, doi:10.1038/nmat2703.

[12] R.J. Harrison, S.A. Redfern, Phys. Earth Planet. Inter. 134 (2002) 253–272, doi:10.1016/S0031-9201(02)00190-5.

[13] R.J. Harrison, S.A.T. Redfern, E.K.H. Salje, Phys. Rev. B 69 (2004) 144101, doi:10.1103/PhysRevB.69.144101.

[14] N. Orlovskaya, Y. Gogotsi, M. Reece, B. Cheng, I. Gibson, Acta Mater. 50 (2002) 715–723, doi:10.1016/S1359-6454(01)00382-2.

[15] N. Orlovskaya, N. Browning, A. Nicholls, Acta Mater. 51 (2003) 5063–5071, doi:10.1016/S1359-6454(03)00354-9.

[16] A.E. Ringwood, Geochim. Cosmochim. Acta 55 (1991) 2083–2110, doi:10.1016/0016-7037(91)90090-R.

[17] T. Yagi, H.-K. Mao, P.M. Bell, Phys. Chem. Miner. 3 (1978) 97–110.

[18] G. Fiquet, A. Dewaele, D. Andrault, M. Kunz, T.L. Bihan, Geophys. Res. Lett. 27 (2000) 21–24.

[19] A.R. Oganov, J.P. Brodholt, G. Price, Earth Planet. Sci. Lett. 184 (2001) 555–560, doi:10.1016/S0012-821X(00)00363-0.

[20] A. Yeganeh-Haeri, D.J. Weidner, E. Ito, Science 243 (4892) (1989) 787–789, doi:10.1126/science.243.4892.787.

[21] S.A.T. Redfern, J. Phys. Condens. Matter 8 (1996) 8267–8275, doi:10.1088/0953-8984/8/43/019.

[22] Z. Wang, S.-I. Karato, K. Fujino, Phys. Earth Planet. Interface 79 (1993) 299–312.

[23] S. Beauchesne, J. Poirier, Phys. Earth Planet. Interface. 55 (1989) 187–199, doi:10.1016/0031-9201(89)90242-2.

[24] R.M. Wentzcovitch, N.L. Ross, G. Price, Phys. Earth Planet. Interface 90 (1995) 101–112, doi:10.1016/0031-9201(94)03001-Y.

[25] M. Alfredsson, J.P. Brodholt, D.P. Dobson, A.R. Oganov, C.R.A. Catlow, S.C. Parker, G.D. Price, Phys. Chem. Miner. 31 (2005) 671–682.

[26] Y. Chen, A. Chernatynskiy, D. Brown, P.K. Schelling, E. Artacho, S.R. Phillpot, Phys. Earth Planet. Interface 210–211 (2012) 75–89, doi:10.1016/j.pepi.2012.08.002.

[27] P. Hirel, Comput. Phys. Commun. 197 (2015) 212–219.

[28] S.J. Plimpton, J. Comp. Phys. 117 (1995) 1–19.

[29] K. Momma, F. Izumi, J. Appl. Crystallogr. 44 (2011) 1272–1276, doi:10.1107/S0021889811038970.

[30] D. Ferré, P. Carrez, P. Cordier, Phys. Earth Planet. Interface 163 (2007) 283–291.

[31] K. Gouriet, P. Carrez, P. Cordier, Model. Simul. Mater. Sci. Eng. 22 (2014) 25020.

[32] H. Jónsson, G. Mills, K.W. Jacobsen, in: Classical and Quantum Dynamics in Condensed Phase Simulations, WORLD SCIENTIFIC, 1998, pp. 385–404, doi:10.1142/9789812839664_0016.

[33] G. Henkelman, H. Jónsson, J. Chem. Phys. 113 (2000) 9978.

[34] A.M. Dziewonski, D.L. Anderson, Phys. Earth Planet. Interface 25 (1981) 297–356, doi:10.1016/0031-9201(81)90046-7.

[35] E. Salje, Annu. Rev. Matter Res. 42 (2012) 1.1–1.19, doi:10.1146/annurev-matsci-070511-155022.

[36] A.R. Lim, J. Phys.: Condens. Matter 7 (1995) 7309, doi:10.1088/0953-8984/7/37/005.