Ultimate mechanical properties of enstatite

Karine Gouriet, Pascal Roussel, Philippe Carrez, Patrick Cordier

► To cite this version:

Karine Gouriet, Pascal Roussel, Philippe Carrez, Patrick Cordier. Ultimate mechanical properties of enstatite. Physics and Chemistry of Minerals, 2022, Physics and chemistry of minerals, 49 (8), pp.30. 10.1007/s00269-022-01206-5. hal-03811834

HAL Id: hal-03811834

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

Submitted on 12 Oct 2022

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
# Ultimate mechanical properties of enstatite

Karine Gouriet<fnref n="1" /> Ⓒ · Pascal Roussel<fnref n="2" /> · Philippe Carrez<fnref n="1" /> · Patrick Cordier<fnref n="1" />,<fnref n="3" />

Received: 16 July 2021 / Accepted: 28 May 2022 / Published online: 11 July 2022
© The Author(s) 2022

## Abstract

The ultimate mechanical properties of MgSiO₃ orthoenstatite (OEN), as characterized here by the ideal strengths, have been calculated under tensile and shear loadings using first-principles calculations. Both ideal tensile strength (ITS) and shear strength (ISS) are computed by applying homogeneous strain increments along high-symmetry directions ([100], [010], and [001]) and low index shear planes ((100), (010), and (001)) of the orthorhombic lattice. We show that the ultimate mechanical properties of OEN are highly anisotropic during tensile loading, with ITS ranging from 4.5 GPa along [001] to 8.7 GPa along [100], and quite isotropic during the shear loading with ISS ranging from 7.4 to 8.9 GPa. During tensile test along [100] and [001], a modified structure close to OEN has been found. This modified structure is more stable than OEN under stress (or strain). We have characterized its elastic and ultimate properties under tensile loading. With ITS ranging from 7.6 GPa along [010] to 25.6 GPa along [001], this modified structure appears to be very anisotropic with exceptional strength along [001].

**Keywords** Orthoenstatite · Ultimate mechanical properties · Ideal tensile strength · Ideal shear strength

## Introduction

The mechanical properties of crystals are usually controlled by defects. In the brittle field microcracks, pores, inclusions control toughness. In the ductile field, flow results from the motion of vacancies, dislocations, twins, disclinations, disconnections, etc. However, in some circumstances, these mechanisms cannot be activated and the mechanical properties of solids are ultimately controlled by the resistance

of bonds. This is the case when plastic deformation mechanisms are inhibited by temperature being too low, or by kinetics reasons at high strain rates. The first approach to this problem was formulated by Frenkel (1926) who predicts the ideal shear stress to be on the order of $G/2\pi$ ($G$ is the shear modulus). Ideal (shear or tensile) stresses are difficult to approach experimentally, but they are now easily accessible using pseudopotential total energy calculation based on the density functional theory (DFT) (Sob et al. 1997; Roundy et al. 1999; Ogata et al. 2004). In simple solids, the correlation between ultimate and elastic properties is quite good, thus validating Frenkel's theoretical approach (Ogata et al. 2004). However, with more complex structures, deviations may occur. Jiang and Srinivasan (2013) have shown that in Fe₃C new bonds are formed during deformation leading to a strong strain-stiffening of the structure. Only calculations at the atomic scale can predict the ultimate properties of these solids.

(Mg, Fe)SiO₃ is one of the most abundant planetary materials, present in planetary mantles, but also in interplanetary dust particles (Rietmeijer, 1998) and chondritic meteorites (Brearley and Jones 1998). During the formation of the solar system, the minerals that make up these objects were subjected to multiple collisions that could lead to the destabilization of the crystalline structure. Hernandez et al. (2020) have recently shown disordering in

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

Karine Gouriet
karine.gouriet@univ-lille.fr

Pascal Roussel
Pascal.Roussel@univ-lille.fr

Philippe Carrez
Philippe.Carrez@univ-lille.fr

Patrick Cordier
patrick.cordier@univ-lille.fr

<footnote n="1">¹ University Lille, CNRS, INRAE, Centrale Lille, UMR 8207-UMET-Unité Matériaux et Transformations, 59000 Lille, France</footnote>

<footnote n="2">² University Lille, CNRS, Centrale Lille, University Artois, UMR 8181, Unité de Catalyse et Chimie du Solide, 59000 Lille, France</footnote>

<footnote n="3">³ Institut Universitaire de France, 1 rue Descartes, 75005 Paris, France</footnote>
shock-compressed enstatite. However, contrary to olivine (Gouriet et al. 2019; Misawa and Shimoyo 2020), the ultimate properties and mechanical stability of enstatite have not yet been documented.

In this study, we investigate the elastic properties and mechanical stability of iron-free MgSiO₃ orthoenstatite. From first-principles calculations, the ideal tensile strengths (ITS) and ideal shear strengths (ISS) are computed along high-symmetry directions [100], [010], and [001] and for homogeneous shear on (100), (010), and (001) planes (here given with respect to the *a*~18.4 Å, *b*~8.9 Å, *c*~5.2 Å cell parameters' setting of orthoenstatite).

## Methods

We performed the derivation of the anisotropic ideal strengths of OEN using the ADAIS free software, written by Zhang et al. (2019). The ADAIS software allows us to implement a homogenous deformation to standard first-principles VASP calculations (Kress and Frurthmüller 1996). Following the extended work of Li et al. (2014) on enstatite polymorphs, all calculations are based on DFT using a plane-wave basis set and the projector augmented wave method (PAW) (Perdew and Wang 1992). The Perdew-Wang (PW91) gradient-corrected functional (GGA) (Wang and Perdew 1991) is employed to take into account the exchange-correlation energy. To ensure suitable atomic force convergence, we used a kinetic energy cutoff of 520 eV for the expansion of the plane-wave basis set. In this work, we use a single grid of 2×4×6 for the k-points sampling, according to a Monkhorst and Pack scheme (Monkhorst and Pack 1976) corresponding to 18 k-points per enstatite (MgSiO₃) unit cell.

Starting from a fully optimized unit cell of MgSiO₃ OEN, we performed tensile and shear tests by applying a homogeneous strain incrementally by step of 0.5%. In practice, crystal atomic layers are moved along the tensile or shear directions (Fig. 1). At each deformation stage, the cell shape and the atomic positions are relaxed, until all the components of the stress tensor are brought to zero, except for the one corresponding to the applied stress condition.

## Results and discussion

### Ground properties

Prior to mechanical testing, the unit cell of OEN (containing 80 atoms) is optimized. The ground-state properties, equilibrium lattice parameters, and elastic constants are given in Tables 1 and 2, respectively. Lattice parameters are found in agreement with experimental values (Jackson et al.

<figure><img src="image_3.png" /><figcaption>**Fig. 1** Illustration of the various loading conditions applied in this study: (a) tensile deformation along [001] where $\epsilon_{zz} = \epsilon$; (b) simple shear [001](100) where $\epsilon_{xz} = \epsilon/2$ and $\epsilon_{zx} = 0$. The image in gray is the structure at $\epsilon=0$, in color at a given $\epsilon$. Mg is in orange, Si in blue, and O in red</figcaption></figure>

2007; Weidner and Vaughan 1982). The elastic constants calculated here for OEN are consistent with those published by Li et al. (2014). Nevertheless, one may notice the well-documented GGA approximation tendency to underestimate the elastic constants with respect to their experimental counterparts (Table 2). As discussed later, it is worth to recall that one of the particularities of the OEN structure consists in alternating SiO₄ tetrahedra chains connected to octahedral Mg1 and Mg2 sites layers along the c-axis and in (100) layers (Periotto et al. 2012) (see suppl. Figure 1). Moreover, one notes that the octahedron on Mg2 site is larger and distorted, compared to the Mg1 octahedron (see suppl. Figure 1b).

### Ideal tensile test calculations

Tensile tests are performed along the [100], [010], and [001] directions. The evolution of the total strain energy as a function of the engineering strain is shown in Fig. 2. First, we observe for all curves a parabolic evolution, which corresponds to the storage of the elastic energy. These parabolic parts correspond to the linear portion on the stress-strain curves shown in Fig. 3.

For both [100] and [001] tensile tests, we observe a discontinuity in the energy curves near 5% strain. However, for the test along [001], after this first discontinuity, the energy curve increases again until it bends at energy values much higher than those reached in the other two directions. The evolution of the energy for the [100] solicitation presents the same feature followed by a monotonous increase until a second energy drop, but of a larger amplitude. The test along [010] shows only one energy drop at 12%, after which the energy decreases continuously with strain, suggesting that above 12% strain, the OEN structure is no longer stable.
<table><caption>**Table 1 Orthoenstatite**</caption><thead><tr><th></th><th>a (Å)</th><th>b (Å)</th><th>c (Å)</th><th>V (Å<sup>3</sup>)</th></tr></thead><tbody><tr><td>This study</td><td>18.399</td><td>8.919</td><td>5.235</td><td>859.05</td></tr><tr><td>Calculated GGA (Li et al. 2014)</td><td>18.417</td><td>8.9345</td><td>5.2448</td><td>863.038</td></tr><tr><td>Calculated LDA (Li et al. 2014)</td><td>18.072</td><td>8.6889</td><td>5.1194</td><td>803.897</td></tr><tr><td>Experimental (Periotto et al. 2012)</td><td>18.2104</td><td>8.820</td><td>5.1767</td><td>821.4</td></tr><tr><td>Experimental (Duffy and Vaughan, 1988)</td><td>18.204</td><td>8.814</td><td>5.176</td><td>830.4</td></tr></tbody></table>

Crystallographic data for OEN at ambient conditions (RT and 0 GPa) compared with values calculated at 0 K and 0 GPa

<table><caption>**Table 2 Orthoenstatite**</caption><thead><tr><th></th><th>C<sub>11</sub></th><th>C<sub>22</sub></th><th>C<sub>33</sub></th><th>C<sub>44</sub></th><th>C<sub>55</sub></th><th>C<sub>66</sub></th><th>C<sub>12</sub></th><th>C<sub>13</sub></th><th>C<sub>23</sub></th></tr></thead><tbody><tr><td>Jackson et al. (2007)</td><td>225</td><td>178</td><td>214</td><td>78</td><td>76</td><td>82</td><td>72</td><td>54</td><td>53</td></tr><tr><td>Weidner and Vaughan (1982)</td><td>236</td><td>173</td><td>216</td><td>84</td><td>79</td><td>80</td><td>74</td><td>57</td><td>50</td></tr><tr><td>Calculated GGA (Li et al. 2014)</td><td>200</td><td>146</td><td>147</td><td>76</td><td>63</td><td>70</td><td>59</td><td>38</td><td>29</td></tr><tr><td>Calculated LDA (Li et al. 2014)</td><td>251</td><td>195</td><td>241</td><td>89</td><td>80</td><td>86</td><td>90</td><td>77</td><td>71</td></tr><tr><td>Orthoenstatite GGA (this study)</td><td>205</td><td>147</td><td>116</td><td>70</td><td>65</td><td>74</td><td>60</td><td>37</td><td>27</td></tr></tbody></table>

Elastic constants (in GPa) for OEN at 0 GPa compared with previous works

The Cauchy stress can be obtained from the derivation of the energy (solid line in Fig. 3). It agrees with the stress state (according to the Hellman–Feynman theorem) of the strained volume (symbols in Fig. 3). The maximums of the stress–strain curves correspond to an inflexion point of the energy–strain curves. Unsurprisingly, the stress–strain curves (Fig. 3) highlight the discontinuities observed on the energy curves. If generally, the maximum of the stress–strain curves corresponds to the ideal strength, it is necessary to understand the origin of these incidents before discussing the observed mechanical properties. To do this, we study the evolution of the structure during deformation, taking the Mg–O bonds as an indicator.

## Evolution of bonding and structural changes

The tensile test along [010] exhibits the simplest behavior. During [010] tensile loading, only the four longest Mg–O bonds are affected: two for the Mg1 site (Fig. 4a/c) and two for the Mg2 site (Fig. 4b/d), see Suppl. Figure 1 for the definition of each site in the enstatite structure. The bond length evolution is small until the maximum value of the stress is reached, where the four Mg–O-bond lengths diverge, as the instability of the structure appears. All equivalent Mg1 and Mg2 exhibit the same behavior. The Mg2 site corresponds to the most distorted octahedron. Then, the longest Mg–O bond (here the bond between Mg2 and O3B) breaks first. Pulling enstatite along [010] leads to a simple behavior with an elastic loading until mechanical instability by Mg–O-bond breaking. The corresponding ITS is 7.6 GPa at 11% strain.

For the tensile test along the [001] direction, the behavior of Mg in site 2 sheds some light on the event observed

at 5% strain where we observe the complete loss of one bond simultaneously with another being re-built (Fig. 5b). Here, bonds of two O3 type oxygen atoms in the B layer are exchanged due to the tilting of the tetrahedra (Figs. 5d, 6). The deformation is accommodated by the elongation of the chain resulting from this tilting (Fig. 6b). This modified structure becomes progressively unstable as, at 16% strain, two Mg1–O-bond lengths increase rapidly and one breaks after 20% strain (Fig. 5a). Then, at 26.5%, an Mg1 loses the second bond, and finally, at 32.5%, an octahedron of the Mg2 site loses one oxygen bond.

The analysis of the tensile test along the [100] direction leads to similar conclusions, with after 8% of strain an exchange of oxygen for the Mg in the site 2 (Fig. 7b) and the formation of the same modified structure as the one described in Fig. 6. After this bond switch, the Mg–O bonds are maintained until the instability occurs at 12%, where four Mg–O bonds are broken (Fig. 7). Here, the distorted octahedron plays an important role in the accommodation of the deformation. After 12%, the modified enstatite structure is finally destabilized.

In Fig. 8, we plot the ratio of the longitudinal over the transversal strains and their evolutions with the engineering strain. We see that this parameter, which is an apparent Poisson ratio, although it is used here beyond its natural field of application which is the linear elasticity, represents a good indicator of structural changes. For example, a clear jump is observed with the exchange of O bonds for Mg2 (Fig. 8a, c) with a marked evolution preceding the structure change (especially for the tensile test along [001]). Comparatively, the evolution is less pronounced at the approach of instability. We also observe several occurrences of the apparent
<figure><img src="image_4.png" /><figcaption>**Fig. 2** Evolution of the total energy as a function of the engineering strain. The tensile directions are (a) [100] (black empty square), (b) [010] (red empty circle), and (c) [001] (blue empty triangle). Insert plots shown a zoom of the beginning of [001] and [001] curves</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Fig. 3** Evolution of the stress as a function of the engineering strain. The tensile directions are [100] (black empty square), [001] (blue empty triangle), and [010] (red empty circle). The instability occurs at the first incident on the curves</figcaption></figure>

Poisson ratio becoming negative, which indicates that enstatite becomes rapidly auxetic under anisotropic strain.

**Modified enstatite: structure and ultimate mechanical properties**

We have further investigated the modified structure obtained at 5–6% strain during tension along [001] and [100]. Once formed, it can be kept when stress is released, as a metastable structure. Indeed, as expected, without applied strain, its energy is larger ($ΔE=0.121$ eV/unit cell) than the one of OEN which is the stable polymorph of $MgSiO_3$ under these conditions. However, above 5% tensile strain along [001] and 8% strain along [100], this modified structure is more stable as shown by the energy drop (Fig. 2).

Even if, as usual, the modified structure has been relaxed without symmetry constraint in the cell, an analysis of the

<figure><img src="image_6.png" /></figure>
<figure><img src="image_7.png" /><figcaption>**Fig. 4** Typical Mg–O-bond length evolution as a function of strain for tensile deformation of the enstatite along [010] for Mg1 (a) and Mg2 (b) sites. To help the reader, we plot also the stress–strain curves. (c)</figcaption></figure>

Mg1 (d) Mg2 octahedra showing the oxygen label and the initial Mg–O-bond lengths

"missed symmetry elements" of the resulting relaxation was undertaken using the ADDSYM tool, as implemented in the Platon software (Spek 2020), a significantly extended/modified implementation of the MISSYM algorithm (Le Page 1987) for the detection of possibly MISsed ADDitional SYMmetry in a given coordinate set. It turns out that this modified structure also adopts the *Pbca* space group of OEN, with relaxed lattice parameters: *a* = 18.727 Å, *b* = 8.837 Å and *c* = 5.363 Å. The modified structure retains the main particularity of the OEN structure, consisting in alternating SiO₄ tetrahedron chains and octahedral Mg1 and Mg2 sites' layers along the *c*-axis and in (100) layers. However, the octahedron on Mg2 site of the modified structure is less distorted than the one of OEN, since, due to the exchange of two O3B oxygens, the bond lengths of both new O3B and O3A are equivalent and close to 2.35 (Figs. 4d and 5d and Table 3), compared to 2.54 and 2.31, respectively. For the octahedron on Mg1 site, almost no change is observed, except the lengths of the two longest bonds that increase by 5% after formation of the modified structure (Figs. 4c

and 5c and Table 3). Concerning the silicon surrounding, if on one hand, the tetrahedral chain A does not suffer from the deformation (the characteristic tetrahedron angle (O3A–O3A–O3A) (Cameron and Papike 1981) changes only by 2.5°); on the other hand, a strong modification is observed for the tetrahedron chain B, due to the lengthening of the chain. The O3B–O3B–O3B angle goes from 139.9° to 158.2° (Supp. Figure 3). In summary, after deformation, both tetrahedron chains tend towards the same configuration, with characteristic angles for the tetrahedral chain close to 160°. Finally, calculated powder X-ray diffractions using Cu Kα₁ radiation have been calculated for OEN and the modified structure, and are given in Supp. Figure 4. It is interesting to note a strong difference between the two diffractograms, giving hope to strong changes in the diffraction pattern, and hence an easy signature to evidence a transition, if an experimental strain is applied along this direction. The modified structure is different from the high-pressure structures found by Jahn (2008) with first principle calculations, called HP-OEN1 in Jahn (2008) or *Pbca*-II in Li et al.
<figure><img src="image_8.png" /><figcaption>**Fig. 5** Typical Mg–O-bond length evolution as a function of strain for tensile deformation along [001] for Mg1 (a) and Mg2 (b) sites. To help the reader, we plot also the stress–strain curves. The gray zone</figcaption></figure>

corresponds to the modified structure. (c) Mg1 and (d) Mg2 octahedra showing the oxygen label and the initial Mg–O-bond lengths, for the relaxed modified structure at 0% of strain

<figure><img src="image_9.png" /><figcaption>**Fig. 6** Chain B of tetrahedral sites: (a) without deformation, enstatite structure (b) modified structure at 6% of strain for tensile deformation along [001]. The Mg–O bonds in green correspond to the lost bond after 6% of deformation. The Mg–O bonds in purple correspond to the new Mg–O bond. The Mg are in orange, O in red, and Si in blue</figcaption></figure>

(2014). Furthermore, at 0 GPa, the difference of internal energy per unit cell between the HP-OEN1 structure found in Jahn (2008) and OEN is below 0.1 eV/unit cell, whereas, here, we computed for the modified structure an increase of the internal energy of 0.121 eV/unit cell. If the structural modifications appear in layer B in the modified structure, only the layer A is affected in HP-OEN1 (Jahn 2008).

Since we can preserve this modified structure in a metastable state, and relax it at 0% strain, we can study its elastic properties. The nine elastic constants computed for the modified structure are given in Table 4. It is worth noticing that all $C_{ij}$ are positive and meets the criteria for mechanical stability of orthorhombic structure, confirming, therefore, the metastability of the MS phase. The results show a larger value of $C_{33}$ than OEN and more isotropic shear moduli (Table 4). We ran also tensile tests along the three directions with this modified structure after relaxation (Fig. 9) which reproduce very well the behavior obtained previously after 4–5% strain (Supp. Figure 2), demonstrating that this portion

<figure><img src="image_10.png" /></figure>
<figure><img src="image_11.png" /><figcaption>**Fig. 7** Typical Mg–O-bond length evolution as a function of strain for tensile deformation along [100] for Mg1 (a) and Mg2 (b) sites. To help the reader, we plot also the stress–strain curves. The gray zone corresponds to the modified stress–strain curve</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Fig. 8** Apparent Poisson ratio (line) as a function of the engineering strain calculated using a unit cell of forsterite under tensile deformation for (a) along [100], (b) along [010], and (c) along [001]. To help the reader, we plot also the stress–strain curves (symbols)</figcaption></figure>
<table><caption>**Table 3** Bond length (in Å) of Mg2–O and Mg1–O, and the O3–O3–O3 (°) angles for orthoenstatite (OEN) and the modified structure (MS) at 0 GPa, compared with experimental data for the OEN</caption><thead><tr><th></th><th>OEN</th><th>Periotto et al.<br>(2012)</th><th>MS</th></tr></thead><tbody><tr><td>Mg1-O1A</td><td>2.046</td><td>2.020</td><td>2.070</td></tr><tr><td>Mg1-O1A</td><td><strong>2.185</strong></td><td><strong>2.146</strong></td><td><strong>2.235</strong></td></tr><tr><td>Mg1-O1B</td><td><strong>2.197</strong></td><td><strong>2.162</strong></td><td><strong>2.235</strong></td></tr><tr><td>Mg1-O1B</td><td>2.076</td><td>2.064</td><td>2.076</td></tr><tr><td>Mg1-O2A</td><td>2.023</td><td>2.006</td><td>2.019</td></tr><tr><td>Mg1-O2B</td><td>2.063</td><td>2.049</td><td>2.025</td></tr><tr><td>Mg2-O1A</td><td>2.088</td><td>2.096</td><td>2.063</td></tr><tr><td>Mg2-O1B</td><td>2.060</td><td>2.057</td><td>2.067</td></tr><tr><td>Mg2-O2A</td><td>2.034</td><td>2.032</td><td>2.044</td></tr><tr><td>Mg2-O2B</td><td>2.000</td><td>1.986</td><td>2.049</td></tr><tr><td>Mg2-O3A</td><td><strong>2.313</strong></td><td><strong>2.297</strong></td><td><strong>2.356</strong></td></tr><tr><td>Mg2-O3B</td><td><strong>2.540</strong></td><td><strong>2.450</strong></td><td><strong>2.343</strong></td></tr><tr><td>O3A–O3A–O3A</td><td>157.6</td><td>159.3</td><td>160.1</td></tr><tr><td>O3B–O3B–O3B</td><td>139.9</td><td>139.4</td><td>158.2</td></tr></tbody></table>

The longest bonds are in bold

of the curves after the incidents describes the behavior of this modified structure and not the one of OEN.

The modified enstatite is very anisotropic with respect to the ideal tensile test investigated. For the [010] tensile test, a small leap appears at 8%, without incident on the structure. The ITS along [100] and [001] are 14.8 and 25.6 GPa, respectively. Along [001], we observe a plateau after 15% strain corresponding to a slight hardening before destabilization of the structure at 30.5%.

The [100] test shows the typical behavior for the stress–strain curve, with a drop just after the maximum of stress. This drop expresses in fact a loss of the modified enstatite structure. Normalizing the ITS by the Young’s modulus shows that the different behaviors between the [100] and [010] directions are mostly due to elasticity, since the normalized values, 0.08 and 0.07, are rather similar (Table 5). However, direction [001] stands out with a normalized ISS about twice as large.

## Ultimate tensile properties of OEN

### Tension

While our tensile test along [010] shows that the mechanical features can be solely attributed to OEN initial structure,

<figure><img src="image_13.png" /><figcaption>**Fig. 9** Modified structure. Stress as a function of the engineering strain calculated using a unit cell of modified enstatite for the tensile tests along [100] direction in black square, [010] direction in red circle, and [001] blue triangles</figcaption></figure>

<table><caption>**Table 5** Modified structure</caption><thead><tr><th>Tensile tests</th><th>[100]</th><th>[010]</th><th>[001]</th></tr></thead><tbody><tr><td>ITS (GPa)</td><td>14.8</td><td>7.6</td><td>25.6</td></tr><tr><td>Corresponding strain (%)</td><td>10.5</td><td>11.0</td><td>23.5</td></tr><tr><td>Young's modulus (GPa)</td><td>176</td><td>102</td><td>186</td></tr><tr><td>Normalized stress</td><td>0.08</td><td>0.07</td><td>0.14</td></tr></tbody></table>

Ideal stresses (and associated engineering strains) determined in this study under tensile for the modified enstatite structure. For tensile tests, we report also the Young’s modulus. The normalized stresses are the ideal stresses divided by the elastic modulus

we have seen that the tests along [001] and [100] quickly involve another structure. The scope of the study must therefore be restricted to the strains which precede the bond switching responsible for the change in structure to describe the mechanical properties of enstatite. The corresponding stress–strain curves are shown in Fig. 10.

Loading enstatite along [001] leads to the smallest ITS, i.e., 4.5 GPa at a critical strain of 5%. This low value is due to the transition to the modified structure which is more stable under strain. The highest stress is sustained when enstatite is loaded along [100], the instability is reached at 6.0% strain for an ITS of 8.7 GPa, here again due to the transition to the modified structure. In between is the ITS value corresponding to [010] tensile loading: 4.5 GPa at 5.0%.

Between the extreme values of the ITS, we have a ratio of 1.9 showing some anisotropy of orthoenstatite consistent with Young’s modulus. For a more consistent comparison, we computed the normalized ITS (ratio between the

<table><caption>**Table 4** Modified structure</caption><thead><tr><th></th><th>C<sub>11</sub></th><th>C<sub>22</sub></th><th>C<sub>33</sub></th><th>C<sub>44</sub></th><th>C<sub>55</sub></th><th>C<sub>66</sub></th><th>C<sub>12</sub></th><th>C<sub>13</sub></th><th>C<sub>23</sub></th></tr></thead><tbody><tr><td>Modified orthoenstatite</td><td>220</td><td>148</td><td>211</td><td>51</td><td>56</td><td>71</td><td>69</td><td>53</td><td>60</td></tr></tbody></table>

Elastic constants (in GPa) for relaxed modified structure at 0 K and 0 GPa at 0% of strain
<figure><img src="image_14.png" /><figcaption>**Fig. 10** Orthoenstatite. Stress as a function of the engineering strain calculated using a unit cell of enstatite for the tensile test along [100] direction in black square, [010] direction in red circle, and [001] blue triangles</figcaption></figure>

ITS and the corresponding Young's modulus). Normalized ITS are comparable (0.05, 0.07, and 0.04, respectively) for [100], [010], and [001] (Table 6). Although pulling along the largest lattice parameter, i.e., along [100], gives the highest Young's modulus (Table 6), the different values for the ITS are only due to elastic behavior. Transverse directions are completely relaxed during the tensile tests, so we can determine Poisson ratios (Table 7) according to the transverse lattice parameters' variations.

## Shear

To complete the ideal properties of OEN, we also performed simple shear deformation tests. Six additional simple shear deformation tests have been performed in this study. The evolutions of stress as a function of the engineering strain are shown in Fig. 11.

In contrast to the previous calculations, more specifically the tensile tests along [100] and [001], we have not seen any formation of the modified structure during the shear tests.

The initial slope of the stress strain curves is consistent with the elastic properties of OEN and we observe that the enstatite structure is maintained until a critical strain, associated with generally the maximum stress value. Ideal shear

<table><caption>**Table 6** Orthoenstatite</caption><thead><tr><th>Tensile tests</th><th>[100]</th><th>[010]</th><th>[001]</th></tr></thead><tbody><tr><td>ITS (GPa)</td><td>8.7</td><td>7.6</td><td>4.5</td></tr><tr><td>Corresponding strain (%)</td><td>6.0</td><td>9.5</td><td>5.0</td></tr><tr><td>Young's modulus (GPa)</td><td>169</td><td>116</td><td>122</td></tr><tr><td>Normalized stress</td><td>0.05</td><td>0.07</td><td>0.04</td></tr></tbody></table>

Ideal stresses (and associated engineering strains) determined in this study under tensile loading for OEN. For tensile tests, we report also the Young's modulus. The normalized stresses are the ideal stresses divided by the Young's modulus

<table><caption>**Table 7** Orthoenstatite</caption><thead><tr><th>Tensile tests</th><th>[100]</th><th>[010]</th><th>[001]</th></tr></thead><tbody><tr><td>ν<sub>[100]</sub></td><td>—</td><td>0.33</td><td>0.14</td></tr><tr><td>ν<sub>[010]</sub></td><td>0.24</td><td>—</td><td>0.09</td></tr><tr><td>ν<sub>[001]</sub></td><td>0.12</td><td>0.09</td><td>—</td></tr></tbody></table>

Poisson ratio determined in this study under tensile tests for strains below 5%

stresses (ISS) range from 7.39 GPa to 8.93 GPa (Table 8). The normalized ISS values (from 0.10 to 0.14) are consistent with ISS value reported for oxides (Ogata et al. 2004) and normalized ISS values obtained for forsterite (Gouriet et al. 2019).

Overall, the stability of the structure mainly depends on Mg–O bonds. Indeed, only the stress–strain curve corresponding to the simple shear test along [100](001) shows a small discontinuity at 6% strain (Fig. 11b). However, the Mg–O-bond analysis of this particular shear test does not show a structural modification. The drop on the stress–strain curve at 6% comes from the Mg2–O3B getting closer to each other. Then, for this [100](001) shear test, a major instability occurs at 20.5% strain when the 2 Mg–O bonds are broken (one bond between Mg1 and O1B (Fig. 12a) and the other between Mg2 and O3B (Fig. 12b)). Consequently, the ISS value of [100](001) shear test reaches 8.93 GPa.

For the other shear tests, two behaviors are observed: one with a small hardening (Fig. 11c) behavior before a jump ([001](010) and [010](001)) and one with only an abrupt jump ([001](100), [100](010) and [010](100)) (see Fig. 11a, b). The Mg–O-bond evolutions are showed in Supp. Figure 5 and 6 for [001](100) and [010](001) respectively.

The jumps after 10% of strain on each stress–strain curves (Fig. 11) are due to the breaking of Mg–O bonds: one bond of each Mg in site 1 (Mg1–O1: Suppl. Figure 5a and Suppl. Figure 6a) and two bonds for Mg in site 2 (Mg2–O3A and O3B: Suppl. Figure 5b and Suppl. Figure 6b). When the Mg in site 2 and the both O3 (each in A and B layers) are separated, we observe the O3B exchange with another one. However, at this stage, the OEN is already lost. The hardening behavior appears when the Mg2–O3B becomes shorter and then more stable.

## Concluding remarks

In this work, the main objective is to investigate the mechanical response of OEN to applied strains until it becomes mechanically unstable. This point is reached after the linear elastic regime, which is exhibited by a parabolic regime in all energy–strain curves. In consequences, in all stress–strain curves, a linear regime is observed. This allows us to
<figure><img src="image_15.png" /><figcaption>**Fig. 11** Orthoenstatite. Stress evolution as a function of the engineering strain calculated using a unit cell of enstatite under shear deformation (a) for simple shear along [100](010) and [010](100), (b) for</figcaption></figure>

simple shear along [100](001) and [001](100), and (c) for simple shear along [001](010) and [010](001) which is shifted to help the reader

<table><caption>**Table 8** Orthoenstatite</caption><thead><tr><th>Simple shear tests</th><th>[010](001)</th><th>[001](010)</th><th>[100](001)</th><th>[001](100)</th><th>[010](100)</th><th>[100](010)</th></tr></thead><tbody><tr><td>ISS (GPa)</td><td>7.48</td><td>7.70</td><td>8.93</td><td>7.39</td><td>7.41</td><td>8.79</td></tr><tr><td>Corresponding strain (%)</td><td>16</td><td>15.0</td><td>20.5</td><td>13.5</td><td>13.5</td><td>16.0</td></tr><tr><td>Shear modulus (GPa)</td><td>73.3</td><td>73.2</td><td>66.0</td><td>68.5</td><td>68.8</td><td>68.9</td></tr><tr><td>Normalized stress</td><td>0.11</td><td>0.10</td><td>0.14</td><td>0.11</td><td>0.13</td><td>0.12</td></tr></tbody></table>

Ideal stresses (and associated engineering strains) determined in this study under simple shear loading. The normalized stresses are the ideal stresses divided by the shear modulus

determine elastic moduli which, with the calculated lattice parameters, validate our calculations.

OEN is quite isotropic under shear loading with ISS value normalized by the shear modulus close to 0.12. The structural analysis shows that the instability occurs when two Mg1–O and Mg2–O bonds are broken, in particular the bond of one Mg with an O3.

During tensile tests, the first stress–strain curve maximums correspond to the ITS of OEN, which are 8.7, 7.6, and 4.5 GPa along the [100], [010], and [001] directions, respectively. The ITS normalized by the Young’s modulus are close to 0.05. During tensile loading along the shortest and the longest lattice parameters, instability appears rapidly during the test (close to 5% of strain). It is important
<figure><img src="image_16.png" /><figcaption>**Fig. 12** Typical Mg–O-bond length evolution as a function of strain for shear deformation along [100](001) for Mg1 (a) and Mg2 (b) sites. To help the reader, we plot also the stress–strain curves</figcaption></figure>

to note that the instability occurs when the system undergoes structural modifications, like broken bonds, leading to structural reconstruction (Jiang and Srinivasan 2013).

Analysis of bond and structure evolution shows that under tensile loading along [100] and [001], OEN gives rise to a modified structure. Under strain, the tetrahedron chain B is stretched, and the O3B–O3B–O3B angle become equivalent to the O3A–O3A–O3A angle in chain A (~160°). Consequently, two O3B are exchanged with the Mg in site 2. This oxygen-bond exchange leads to the modified structure. We have characterized some properties of the modified structure by computing the elastic constants, ITS, and the powder X-ray diffraction pattern.

This modified structure becomes more resistant against deformation, with two larger ITS values. Along the [100] direction, the strain reaches 10.5% before the instability occurs, with the ITS value being 1.7 times larger than the ITS value of OEN. Along the [001] direction, the instability is reached at 23.5% of strain, which corresponds to an ITS value 5.7 times larger than the ITS value of OEN. Along the [010] direction, the two structures exhibit the same ITS value with comparable ultimate strains (9.5–11%).

After the case of cementite described by Jiang and Srinivasan (2013), this study shows that in OEN, large strains can lead to Mg–O-bond rearrangements responsible for a significant strain-stiffening. In cementite, this behavior was observed under shear loading. Here, we report the first occurrence of strain-stiffening under tensile loading and we show that this behavior is due to the formation of a modified structure which can be brought at ambient conditions (no applied stress) in a metastable state. Hence, some physical properties of this modified structure are presented.

**Supplementary Information** The online version contains supplementary material available at https://doi.org/10.1007/s00269-022-01206-5.

**Author contributions** Conceptualization, PC; methodology, PC; validation, PC; investigation, KG; data analysis, KG and PR; writing—original draft preparation, KG; writing—review and editing, PR, PC, and PC; funding acquisition, PC.

**Funding** This research was funded by European Research Council (ERC) under the European Union’s Horizon 2020 Research and Innovation Program under grant agreement No 787198—TimeMan.

**Availability of data and materials** CIF-files of OEN and the modified structure can be uploaded.

**Code availability** Not applicable.

## Declarations

**Conflict of interest** Not applicable.

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
## References

Brearley AJ, Jones RH (1998) Chondritic meteorites. In: Papike JJ (ed) Planetary materials reviews in mineralogy, vol 36. Walter de Gruyter GmbH & Co KG, Berlin, p C1

Cameron M, Papike JJ (1981) Structural and chemical variations in pyroxenes. Am Miner 66:1–50

Duffy TS, Vaughan MT (1988) Elasticity of enstatite and its relationship to crystal structure. J Geophys Res Solid Earth 93(B1):383–391

Frenkel J (1926) Zur Theorie der Elastizitätsgrenze und der Festigkeit kristallinischer Körper. Z Phys 37:572–609

Gouriet K, Carrez P, Cordier P (2019) Ultimate mechanical properties of forsterite. Minerals 9(12):787

Hernandez J-A, Morard G, Guarguaglini M, Alonso-Mori R, Benuzzi-Mounaix A, Bolis R et al (2020) Direct observation of shock-induced disordering of enstatite below the melting temperature. Geophys Res Lett 47:e2020GL088887

Jackson J, Sinogeikin SV, Bass JD (2007) Sound velocities and single-crystal elasticity of orthoenstatite to 1073 K at ambient pressure. Phys Earth Planet Int 161:1–12

Jahn S (2008) High-pressure phase transitions in MgSiO₃ orthoenstatite studied by atomistic computer simulations. Am Miner 93:528–532

Jiang C, Srinivasan SG (2013) Unexpected strain-stiffening in crystal-line solids. Nature 496:339–342

Kresse G, Furthmuller J (1996) Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set. Phys Rev B 54(16):11169–11186

Le Page Y (1987) Computer derivation of the symmetry elements implied in a structure description. J Appl Crystallogr 20:264–269

Li B, Kung J, Wei L, Liebermann R (2014) Phase transition and elasticity of enstatite under pressure from experiments and first-principles studies. Phys Earth Planet Int 228:63–74

Misawa M, Shimoyo F (2020) First-principles study of pressure-induced amorphization of Fe₂SiO₄ fayalite. Phys Status Solidi B 257:2000173

Monkhorst HJ, Pack JD (1976) Special points for Brillouin-zone integrations. Phys Rev B 13(12):5188–5192

Ogata S, Li J, Hirosaki N, Shibutani Y, Yip S (2004) Ideal shear strain of metals and ceramics. Phys Rev B 70:104404

Perdew JP, Wang Y (1992) Accurate and simple analytic representation of the electron-gas correlation-energy. Phys Rev B 45(23):13244–13249

Periotto B, Balic-Zunic T, Nestola F, Katerinopoulou A, Angel R (2012) Re-investigation of the crystal structure of enstatite under high-pressure conditions. Am Miner 97:1741–1748

Rietmeijer FJM (1998) Interplanetary dust particles. In: Papike JJ (ed) Planetary materials reviews in mineralogy, vol 36. Walter de Gruyter GmbH & Co KG, Berlin, p B1

Roundy D, Krenn CR, Cohen ML, Morris JW Jr (1999) Ideal shear strengths of fcc aluminum and copper. Phys Rev Lett 82:2713–2716

Spek AL (2020) CheckCIF validation ALERTS: what they mean and how to respond. Acta Crystallogr A E76:1–11

Sob M, Wang LG, Vitek V (1997) Theoretical tensile stress in tungsten single crystals by full-potential first-principles calculations. Mater Sci Eng A234–236:1075–1078

Wang Y, Perdew JP (1991) Correlation hole of the spin-polarized electron gas, with exact small-wave-vector and high-density scaling. Phys Rev B 44(24):13298–13307

Wang J, Zhou Y, Lin ZL (2007) Mechanical properties and atomic deformation mechanism of γ-Y₂Si₂O₇ from first-principles investigations. Acta Mater 55:6019–6026

Weidner DJ, Vaughan MT (1982) Elasticity of pyroxenes: effects of composition versus crystal structure. J Geophys Res Solid Earth 87(B11):9349–9353

Zhang SH, Fu ZH, Zhang RF (2019) ADAIS: Automatic derivation of anisotropic ideal strength via high-throughput first-principles computations. Comput Phys Commun 238:244–253

**Publisher's Note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.