Ultimate mechanical properties of enstatite

The ultimate mechanical properties of MgSiO3 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]. Supplementary Information The online version contains supplementary material available at 10.1007/s00269-022-01206-5.


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 ( 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 3 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 3 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 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.

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 3 ) unit cell.
Starting from a fully optimized unit cell of MgSiO 3 OEN, we performed tensile and shear tests by applying a homogenous 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.

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. 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 welldocumented 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 4 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. 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.  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.

Evolution of bonding and structural changes
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 Table 1 Orthoenstatite Crystallographic data for OEN at ambient conditions (RT and 0 GPa) compared with values calculated at 0 K and 0 GPa 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 The instability occurs at the first incident on the curves "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 4 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α 1 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 diffractogramms, 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.  (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  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.

Tension
While our tensile test along [010] shows that the mechanical features can be solely attributed to OEN initial structure, 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 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 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)  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     , 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. Availability of data and materials CIF-files of OEN and the modified structure can be uploaded.

Author contributions
Code availability Not applicable.

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:// creat iveco mmons. org/ licen ses/ by/4. 0/.