Nonlinear emergent elasticity and structural transitions of skyrmion crystal under uniaxial distortion

TL;DR

This paper investigates the nonlinear elasticity and structural phase transitions of skyrmion crystals under uniaxial strain, revealing complex behaviors including multiple phase transitions, emergent polymorphism, and metastable states.

Contribution

It provides a theoretical analysis of skyrmion crystal deformation under uniaxial stress using three methods, uncovering new phases and nonlinear elastic properties.

Findings

Skyrmion crystals exhibit negative elastic stiffness under moderate tension.

Up to six phase transitions occur under strong strain, forming four new metastable emergent crystals.

Emergent crystals display global rotation and unlimited deformability, leading to exotic polymorphism.

Abstract

Emergent crystals are periodic alignment of "emergent particles", i.e., localized collective behavior of atoms or their charges/spins/orbits. These novel states of matter, widely observed in various systems, may deform under mechanical forces with elasticity strikingly different from that of the underlying material. However, their nonlinear and critical behaviors under strong fields are hitherto unclear. Here we theoretically study the nonlinear elasticity and structural transitions of skyrmion crystals (SkX) suffering uniaxial distortion by using three different methods. Under moderate tension, SkX behaves like a ductile material, with a negative crossover elastic stiffness and a negative emergent Poisson's ratio at appropriate conditions of magnetic field. Under strong straining, we observe at most six phase transitions, leading to appearance of four novel emergent crystals that are thermodynamically metastable. When subject to external loads, emergent crystals rotate globally, and their composing "particles" have unlimited deformability, which render their exotic polymorphism.

Equations

$$\tilde{w}(\mathbf{m})=$$

$$+t{\mathbf{m}}^2+{\mathbf{m}}^4+\sum _{i=1}^3\left[{\tilde{A}}_c{m}_i^4+{\tilde{A}}_e{\left(\frac{\partial m_i}{\partial r_i}\right)}^2\right]$$

$$+{\epsilon }_{11}[\tilde{K}{\mathbf{m}}^2+{\tilde{L}}_1{m}_1^2+{\tilde{L}}_2{m}_3^2+{\tilde{L}}_{O1}(m_{1,2}{m}_3$$

$$-m_{1,3}{m}_2)+{\tilde{L}}_{O2}(m_{3,1}{m}_2-m_{2,1}{m}_3)$$

$$+{\tilde{L}}_{O3}{m}_1(m_{2,3}-m_{3,2})],$$

$$\mathbf{m}={\mathbf{m}}_0+\sum _{i=1}^n\sum _{j=1}^{n_i}{\mathbf{m}}_{{\mathbf{q}}_{ij}}{e}^{\mathrm{i}{\mathbf{q}}_{ij}({\mathbf{\epsilon }}^{ea})\cdot \mathbf{r}},$$

$${\mathbf{\epsilon }}^{ea}=[{\epsilon }_{11}^e,{\epsilon }_{22}^e,{\epsilon }_{12}^e,{\omega }^e{]}^T,$$

$${\mathbf{m}}^a=[$$

$${\tilde{c}}_{121}^{re},{\tilde{c}}_{122}^{im},{\tilde{c}}_{123}^{re},{\tilde{c}}_{131}^{re},{\tilde{c}}_{132}^{im},{\tilde{c}}_{133}^{re}{]}^T.$$

$$\tilde{w}(\mathbf{m},{\epsilon }_{ij})=$$

$$t{\mathbf{m}}^2+{\mathbf{m}}^4+\sum _{i=1}^3\left[{\tilde{A}}_c{m}_i^4+{\tilde{A}}_e{\left(\frac{\partial m_i}{\partial r_i}\right)}^2\right]$$

$$+{\tilde{w}}_{el}+{\tilde{w}}_{me},$$

$${\tilde{w}}_{el}=$$

$$+{\epsilon }_{22}^2{\epsilon }_{33}^2)+\frac{1}{2}{\tilde{C}}_{44}({\gamma }_{12}^2+{\gamma }_{13}^2+{\gamma }_{23}^2)$$

$${\tilde{w}}_{me}=$$

$$+{\tilde{L}}_2(m_3^2{\epsilon }_{11}+m_1^2{\epsilon }_{22}+m_2^2{\epsilon }_{33})+{\tilde{L}}_3(m_1{m}_2{\gamma }_{12}$$

$$+m_1{m}_3{\gamma }_{13}+m_2{m}_3{\gamma }_{23})+\sum _{i=1}^6{\tilde{L}}_{Oi}{\tilde{f}}_{Oi},$$

$${\tilde{f}}_{O1}=$$

$$-m_{2,1}{m}_3)+{\epsilon }_{33}(m_{3,1}{m}_2-m_{3,2}{m}_1),$$

$${\tilde{f}}_{O2}=$$

$$-m_{3,2}{m}_1)+{\epsilon }_{33}(m_{2,3}{m}_1-m_{1,3}{m}_2),$$

$${\tilde{f}}_{O3}=$$

$$+{\epsilon }_{33}{m}_3(m_{1,2}-m_{2,1}),$$

$${\tilde{f}}_{O4}=$$

$$-m_{2,3}{m}_3)+{\gamma }_{12}(m_{3,2}{m}_2-m_{3,1}{m}_1),$$

$${\tilde{f}}_{O4}=$$

$$-m_{2,3}{m}_3)+{\gamma }_{12}(m_{3,2}{m}_2-m_{3,1}{m}_1),$$

$${\tilde{f}}_{O5}=$$

$$-m_{3,2}{m}_3)+{\gamma }_{12}(m_{2,3}{m}_2-m_{1,3}{m}_1),$$

$${\tilde{f}}_{O6}=$$

$$+{\gamma }_{12}{m}_3(m_{2,2}-m_{1,1}).$$

$$\mathbf{m}={\mathbf{m}}_0+\sum _{i=1}^{\infty }\sum _{j=1}^{n_i}{\mathbf{m}}_{{\mathbf{q}}_{ij}}{e}^{\mathrm{i}{\mathbf{q}}_{ij}\cdot \mathbf{r}}.$$

$${\mathbf{q}}_{11}=$$

$${\mathbf{q}}_{12}=$$

$$-1-{\epsilon }_{11}^e+\sqrt{3}({\epsilon }_{12}^e+{\omega }^e){]}^T,$$

$$\bar{w}=$$

$$+\frac{1}{2}(d\mathbf{\epsilon }{)}^T\tilde{\mathbf{C}}d\mathbf{\epsilon }+(d\mathbf{\epsilon }{)}^T\tilde{\mathbf{h}}d{\mathbf{\epsilon }}^{ea}$$

$$+(d{\mathbf{\epsilon }}^{ea}{)}^T{\tilde{\mathbf{g}}}^{em}d{\mathbf{m}}^a+(d\mathbf{\epsilon }{)}^T{\tilde{\mathbf{g}}}^{m}d{\mathbf{m}}^a,$$

$${\mathbf{\epsilon }}^{ea}=[{\epsilon }_{11}^e,{\epsilon }_{22}^e,{\epsilon }_{12}^e,{\mathbf{\omega }}^e{]}^T,$$

$$\mathbf{\epsilon }=[{\epsilon }_{11},{\epsilon }_{22},{\epsilon }_{12},{\gamma }_{23},{\gamma }_{13},{\gamma }_{12}{]}^T,$$

$${\mathbf{m}}^a=[$$

$${\tilde{c}}_{122}^{im},{\tilde{c}}_{123}^{re},{\tilde{c}}_{131}^{re},{\tilde{c}}_{132}^{im},{\tilde{c}}_{133}^{re},{\tilde{c}}_{211}^{re},\cdot \cdot \cdot {]}^T.$$

$${\tilde{C}}_{ij}={\left(\frac{{\partial }^2\bar{w}}{\partial {\epsilon }_i\partial {\epsilon }_j}\right)}_0,{\tilde{h}}_{ij}={\left(\frac{{\partial }^2\bar{w}}{\partial {\epsilon }_i\partial {\epsilon }_j^{ea}}\right)}_0,$$

$${\tilde{g}}_{ij}^m={\left(\frac{{\partial }^2\bar{w}}{\partial {\epsilon }_i\partial m_j^a}\right)}_0,{\tilde{C}}_{ij}^e={\left(\frac{{\partial }^2\bar{w}}{\partial {\epsilon }_i^{ea}\partial {\epsilon }_j^{ea}}\right)}_0,$$

$${\tilde{g}}_{ij}^{em}={\left(\frac{{\partial }^2\bar{w}}{\partial {\epsilon }_i^{ea}\partial m_j^a}\right)}_0,{\tilde{\mu }}_{ij}^a={\left(\frac{{\partial }^2\bar{w}}{\partial m_i^a\partial m_j^a}\right)}_0.$$

Full text

Introduction

A recent outburst of interest in emergent crystals occurs due to the discovery of skyrmions with various chiralitymuhlbauer2009skyrmion ; kezsmarki2015neel , commensurabilityvon2015influence ; langner2014coupled and dimensionalitymuhlbauer2009skyrmion ; tanigaki2015real , and “relatives” of skyrmionslee2015zero ; huang2017stabilization ; ackerman2015self in different systemsbauerle1996laboratory ; al2001skyrmions ; rossler2006spontaneous ; fu2016persistent ; nych2017spontaneous , especially in bulkmuhlbauer2009skyrmion ; tanigaki2015real ; karube2016robust ; johnson2013mnsb ; leonov2015multiply and geometrically confinedyu2010real ; heinze2011spontaneous ; yu2011near ; jiang2017skyrmions ; du2015electrical magnetic materials where they appear as localized spin textures with nontrivial topology. These topological objects behave as elementary particles in an emerging world, such that they form novel crystalline states at low temperature, and melt under heatingambrose2013melting . On one side, these emergent crystals benefit from the exotic local properties of their composing particles (e.g., topological protection, mobility to electric currentjonietz2010spin ; zang2011dynamics , and skyrmion Hall effectsneubauer2009topological ; litzius2017skyrmion ; jiang2017direct ). On the other side, they give rise to novel macroscopic emergent properties when interacting with external fieldsshiomi2013topological ; white2014electric ; shibata2015large ; okamura2016transition . Systematic study of the deformation and instability of emergent crystals under different types of external fields provides a major challenge for us to understand the fundamental similarities and discrepancies between emergent crystals and ordinary crystals, and to develop reliable approaches toward effective manipulation of these new states of matter.

Skyrmion crystals (SkX) in chiral magnets and their layered structures are profoundly affected by deformation of the material. It is known that uniaxial loadingnii2015uniaxial ; chacon2015uniaxial ; wang2018uniaxial , as well as misfit strains in thin film structureskarhu2010structure ; huang2012extended both significantly change the stability of SkX in the temperature-magnetic field phase diagram. Moreover, mechanical loads affect the generationliu2017chopping and chiralitychen2015unlocking of skyrmion, and the elementary excitationszhang2017ultrasonic in the SkX. A fundamental type of phenomena that connects or even explains the ones mentioned above is the ability of SkX to deform under strains: when subject to uniaxial tension, the SkX inside FeGe thin film is found to undergo significant deformation about 66 times larger than the deformation of the underlying atomic latticeshibata2015large , and this emergent elasticity of SkX has also been observed in MnSishibata2015large ; fobes2017versatile ; kang2017elastic . At present stage, most of existing theories on this anomalous emergent deformation of SkX are restricted to its elastic rangehu2016emergent ; petrova2011spin , while the nonlinear deformation and instability of SkX suffering moderate or strong elastic fields have never been elucidated.

In this work, we study the nonlinear elasticity and structural phase transitions of SkX in bulk helimagnets within an analytical theoretical framework, where the deformation of SkX is described by the emergent elastic strains and the emergent rotation angleshu2016emergent . Developed upon a thermodynamic potential incorporating a comprehensive magnetoelastic functional for B20 helimagnetshu2017unified , the framework has proved its effectiveness by quantitatively reproducing for MnSi the phase diagrammuhlbauer2009skyrmion , the variation of elastic stiffnessnii2014elastic , and the changed stability of SkX in thin filmsli2013robust ; hu2018effect . For SkX in MnSi suffering uniaxial tension, we find that the stress-emergent strain curve before phase transition resembles the stress-strain curve of typical ductile materials (e.g., aluminum alloy). However, in most of the conditions studied, SkX possesses a negative emergent Poisson s ratio. And at high magnetic field and small uniaxial strain, it has an exotic negative crossover elastic stiffness. As the strain varies, the SkX undergoes five or six times of subsequent structural phase transitions, leading to appearance of four distorted SkX phases. Further analysis shows that these new states are thermodynamically metastable. We quantitatively reproduce the exotic emgergent elasticity of SkX observed in FeGe thin film suffering uniaxial tensionshibata2015large , and find that it is caused by a triangle-square phase transition.

Model

Consider a cubic helimagnet stabilized in the SkX phase suffering uniaxial distortion along the x-axis, the equilibrium state of the system is determined by the following rescaled Helmholtz free energy density functional (see the Methods section for derivation)hu2017unified

[TABLE]

where $\mathbf{m}$ is the rescaled magnetization, $\mathbf{\varepsilon}_{11}$ is the strain component describing the uniaxial distortion, $t$ is the rescaled temperature, $\mathbf{b}$ is the rescaled magnetic field, and $\widetilde{w}_{0}$ denotes the part of free energy density independent of the magnetization. To describe the deformed SkX, $\mathbf{m}$ is expanded within the $n^{\text{th}}$ order Fourier representation ashu2016emergent ; hu2018wave

[TABLE]

or alternatively $\mathbf{m}=\mathbf{(}\mathbf{m}^{a},\bm{\varepsilon}^{ea})$, where the vector $\mathbf{m}^{a}$ contains all components of the vectors of Fourier magnitude $\mathbf{m}_{0}$ and $\mathbf{m}_{\mathbf{q}_{ij}}$ for all $i$ and $j$, and the vector $\bm{\varepsilon}^{ea}$ contains all components of emergent elastic strains and the emergent rotation angles. For 2-D SkX. the exact form of $\mathbf{m}^{a}$ and $\bm{\varepsilon}^{ea}$ for $1^{st}$ order Fourier representation can be written as (see the Methods section for derivation)

[TABLE]

[TABLE]

One should keep in mind that the length of $\mathbf{m}^{a}$ depends on the order of Fourier representation, For $2^{nd}$, $3^{rd}$ and $4^{th}$ order Fourier representation, the number of components of $\mathbf{m}^{a}$ becomes 21, 30 and 48, respectively. At given values of $t$, $\mathbf{b}$, and $\varepsilon_{11}$, the equilibrium state of the SkX phase is determined by minimizing $\bar{w}(\mathbf{m}^{a},\bm{\varepsilon}^{ea})$ with respect to all components of $\mathbf{m}^{a}$ and $\bm{\varepsilon}^{ea}$, where $\bar{w}(\mathbf{m}^{a},\bm{\varepsilon}^{ea})=\frac{1}{V}\int_{V}\widetilde{w}(\mathbf{m}^{a},\bm{\varepsilon}^{ea})$ denotes the averaged free energy density of the system. To determine the thermodynamic equilibrium state of the system at given values of $t$, $\mathbf{b}$, and $\varepsilon_{11}$, we have to compare the minimized value $\bar{w}_{min}(\mathbf{m}^{a},\bm{\varepsilon}^{ea})$ of the SkX phase with that of all the other possible magnetic phases such as the ferromagnetic phase and the conical phase. Concerning the large number of variables included in $\mathbf{m}^{a}$ and $\bm{\varepsilon}^{ea}$, the minimization problem is solved numerically. To guarantee that the solution obtained actually corresponds to a local minimum of the free energy of the system, the method of soft-mode analysishu2018wave is incorporated in the minimization process.

Nonlinear emergent elasticity of skyrmion crystal in MnSi

We analyze the nonlinear emergent elasticity and instability of SkX in a prototype helimagnet MnSi (thermodynamic parameters shown in ref.hu2017unified ) by calculating variation of the equilibrium field configuration of SkX (i.e., the values of components of $\bm{\varepsilon}^{ea}$ and $\mathbf{m}^{a}$) with the uniaxial distortion. We assume that the SkX is distributed in the x-y plane and the rescaled magnetic field satisfies $\mathbf{b}=[0,0,b]^{T}$.

We first use three different methods to calculate the emergent elastic property of SkX at small uniaxial tensile strain, including the theory of emergent elasticity, the $4^{th}$ order Fourier representation based free energy minimization, and the Monte Carlo simulation based free energy minimization. As shown by the $\varepsilon_{11}-\varepsilon^{e}_{11}$ curves obtained at $t=0.5$ in FIG. 1(a), the results obtained from the three methods agree quantitatively well in the range $0\leq\varepsilon_{11}\leq 0.0016$, and are exactly equal to each other when $\varepsilon_{11}\to 0$. In this range, $\varepsilon^{e}_{11}$ is almost linear with $\varepsilon_{11}$ at all calculated conditions of $b$ and $\widetilde{A}_{e}$ according to the two free energy minimization methods. In a wider range of $\varepsilon_{11}$ (FIG. 1(a, b)), one discovers that the $\varepsilon_{11}-\varepsilon^{e}_{11}$ curves gradually presents nonlinearity, and the discrepancy between results obtained by the two methods gradually increases while maintaining qualitative agreement. This qualitative agreement repeats itself in the $\varepsilon_{11}-\varepsilon^{e}_{22}$ curves (FIG. 2(a)), demonstrating the effectiveness of the Fourier representation based free energy minimization method.

In the theory of emergent elasticity hu2016emergent , the emergent elastic property of SkX is characterized by the crossover strain ratio (CSR) matrix, whose components are the linear coefficients relating $\bm{\varepsilon}^{ea}$ and the elastic strains of the underlying material. FIG. 1(a, c) shows that the strain-free CSR component $\lambda_{11}$ (inverse of the gradient of the $\varepsilon_{11}-\varepsilon^{e}_{11}$ curves at $\varepsilon_{11}=0$) decreases sharply with $b$, and is slightly affected by a change of exchange anisotropy (from $\widetilde{A}_{e}=0$ to $\widetilde{A}_{e}=-0.05$). Specifically, the $\varepsilon_{11}-\varepsilon^{e}_{11}$ curves calculated at $b=0.3$ correspond to a negative strain-free $\lambda_{11}$, indicating a compressing SkX in the x-axis when the underlying material is suffering elongation in the same direction.

In a wider range $0\leq\varepsilon_{11}\leq 0.02$, the $\varepsilon_{11}-\varepsilon^{e}_{11}$ curves (FIG. 1(b)) exhibit a similar shape to the stress-strain curve of any typical ductile metal (e.g., Aluminum alloy). Multiplying the vertical axis of FIG. 1(b) by the elastic constant $C_{11}$, we obtain the stress-emergent elastic strain curves for SkX under uniaxial distortion. At all calculated conditions of $b$ and $\widetilde{A}_{e}$, the emergent deformation of SkX with $\varepsilon_{11}$ can be divided in the three stages: linear elastic region, nonlinear elastic region, and “failure”, where the strain corresponding to the point separating the linear elastic region and nonlinear elastic region is called emergent yielding strain $\varepsilon_{Y}$, and the strain corresponding to the failure point is called emergent failure strain $\varepsilon_{F}$. Imitating the definition of yielding point for ductile materials such as Aluminum alloy, we define the emergent yielding point by the intersection of the $\varepsilon_{11}-\varepsilon^{e}_{11}$ curve and the straight line determined by the strain-free $\lambda_{11}$ right shifted by $2\%$ on the $\varepsilon^{e}_{11}$-axis. And we define the failure point to be the point where a structural phase transition occurs. $\varepsilon_{Y}$ provides the upper boundary concerning the applicability of the theory of emergent elasticityhu2016emergent , while $\varepsilon_{F}$ denotes the critical strain where a drastic change of properties of SkX is about to occur. From FIG. 1(b, c), we learn that $\varepsilon_{Y}$ decreases with $b$ while $\varepsilon_{F}$ increases with $b$, both of which are induced by the decrease of strain-free $\lambda_{11}$ with $b$. To be more specific, consider the $\varepsilon_{11}-\varepsilon^{e}_{11}$ curve at $b=0.3$, $\varepsilon^{e}_{11}$ first decreases and then increases with $\varepsilon_{11}$, indicating a drastic change of $\lambda_{11}$ with $\varepsilon_{11}$, and hence $\varepsilon_{Y}$ is small since it describes the region where $\lambda_{11}$ approximately keeps its initial strain-free value. Meanwhile, $\varepsilon_{F}$ is large because larger value of $\varepsilon_{11}$ is required to achieve the same level of $\varepsilon^{e}_{11}$ than the case where $b$ is smaller.

FIG. 2 shows the variation of $\varepsilon^{e}_{22}$, $\varepsilon^{e}_{12}$ and $\omega^{e}$ with $\varepsilon_{11}$ in the range $0\leq\varepsilon_{11}\leq 0.02$, where the $\varepsilon_{11}-\varepsilon^{e}_{12}$ curves and $\varepsilon_{11}-\omega^{e}$ curves plotted in FIG. 2(c, d) are calculated at the condition $\widetilde{A}_{e}=-0.05$ (neglecting exchange anisotropy, $\varepsilon^{e}_{12}$ and $\omega^{e}$ are both zero for any $\varepsilon_{11}$). It is shown that $\varepsilon^{e}_{22}$ varies in a complicated way with $\varepsilon_{11}$ (FIG. 2(a)). We define the emergent Poisson s ration as $\nu^{e}=-\varepsilon^{e}_{22}/\varepsilon^{e}_{11}$, and plot in FIG. 2(b) the $\nu^{e}-\varepsilon_{11}$ curves calculated at various conditions. It is found that for all calculated conditions, SkX always possess negative emergent Poisson s ratio in a considerable range of positive $\varepsilon_{11}$. A negative emergent Poisson s ratio indicates the number of skyrmions in a unit volume of the material suffering mechanical forces is not conservative.

Structural transitions of skyrmion crystal in MnSi suffering uniaxial tension

The break points appearing near $\varepsilon_{11}=0.0001$ in the $\varepsilon_{11}-\varepsilon^{e}_{11}$, $\varepsilon_{11}-\varepsilon^{e}_{12}$ and $\varepsilon_{11}-\omega^{e}$ curves at $b=0.1$ shown in FIG. 1, 2 imply the occurrence of a structural phase transition of the SkX, accompanied by emergent rotation (since $\omega^{e}$ changes sharply). By studying the metastability and phase transitions of SkX in MnSi at $t=0.5$, $b=0.3$ suffering uniaxial distortion, we find that $\mathbf{a})$ 6 types of nontrivial emergent crystalline states may appear in MnSi suffering uniaxial tension, distinguished by a difference of emergent elastic property. They are: triangle skyrmion crystal (T-SkX), rotated triangle skyrmion crystal (RT-SkX), deformed square skyrmion crystal (DS-SkX), elongated and rotated triangle skyrmion crystal (ERT-SkX), unstable skyrmion crystal (U-SkX), and in-plane single-Q (IPSQ) phase. The field configuration of these different phases are shown in the insets of FIG. 3. $\mathbf{b}$) For small positive $\varepsilon_{11}$, T-SkX and RT-SkX are the two possible states that may appear, where RT-SkX is obtained by rotating the T-SkX by about $30^{\circ}$. Neglecting exchange anisotropy, T-SkX is the only stable state in this region, which switches to RT-SkX when $\widetilde{A}_{e}=-0.05$. If we change the sign of $\varepsilon_{11}$ to apply a small compression, two structural transitions will swiftly occur in turn, first from T-SkX to RT-SkX (or vice versa) and then to the ferromagnetic phase. $\mathbf{c}$) As $\varepsilon_{11}$ further increases, a triangle-square transition of SkX will always occur leading to appearance of the DS-SkX phase. The local field configuration of the emergent particle in the DS-SkX phase is always notably elongated in the direction of uniaxial distortion (FIG. 3(a, b)). Interestingly, we find that immediately after the phase transition, the change of emergent elastic strain $\Delta\varepsilon^{e}_{11}$ compared with the strain-free SkX is about 60 times larger than the uniaxial strain $\varepsilon_{11}$ applied to the underlying material. This result agrees quantitatively well with the significant deformation of the SkX observed in FeGe thin film suffering uniaxial distortionshibata2015large . Concerning the similarity of crystalline structure between MnSi and FeGe, we argue that this phase transition from RT-SkX to DS-SkX explains the exotic emergent elasticity of SkX in FeGe. $\mathbf{d})$ As $\varepsilon_{11}$ further increases, the DS-SkX will transform to the ERT-SkX phase, where the emergent particle is severely elongated in the direction of uniaxial distortion, and the lattice orientation and direction of elongation is are dominated by the loading direction. Neglecting exchange anisotropy (FIG. 3(c)), the long edges of the elongated skyrmions in the ERT-SkX phase are exactly parallel to the direction of uniaxial distortion, while for $\widetilde{A}_{e}=-0.05$ (FIG. 3(a, b)), these long edges are slightly rotated. We also notice that the location of phase transition point from DS-SkX to ERT-SkX is sensitive to the loading direction (FIG. 3(a, b)). $\mathbf{e)}$ As $\varepsilon_{11}$ further increases, a phase transition from the ERT-SkX phase to the U-SkX phase may occur, depending on the loading direction. When the uniaxial distortion is applied in the x-axis (FIG. 3(a, c)), we see a clear sign of second order phase transition by observing the variation of order parameters with $\varepsilon_{11}$, indicating appearance of the U-SkX phase, which is mechanically very soft and unstable. It collapses into the IPSQ phase easily under further extension. On the other hand, if we apply uniaxial distortion in the y-axis (FIG. 3(b)), ERT-SkX is the only stable state in the same range of $\varepsilon_{11}$. As $\varepsilon_{11}$ increases, the ERT-SkX phase also becomes gradually softer and finally transform to the IPSQ phase. $\mathbf{f)}$ In FIG. 3, we choose $\varepsilon^{e}_{11}$ and $\omega^{e}$ as the order parameters of the system, and plot their variation with $\varepsilon_{11}$ to illustrate the behavior of different phases. A novel behavior of SkX as well as any emergent crystalline states is the general existence and variation of $\omega^{e}$ when subject to mechanical forces. Emergent rotation exists for two major reasons. Firstly, as $\varepsilon_{11}$ increases, the anisotropy of the free energy density functional is gradually changing, which renders a change of lattice direction of SkX to achieve lower free energy density. Secondly, the structures of different emergent crystalline states (e.g., T-SkX and DS-SkX) appearing on the two sides of a transition are usually directionally incompatible with each other. Therefore, emergent rotation is inevitable which leads to jump of $\omega^{e}$ during these phase transitions. Emergent rotation with respect to the underlying atomic lattices is an important and unique way for SkX and other emergent crystals to resist tension of the material. Nevertheless, since emergent rotation breaks the ordinary periodic boundary condition widely used in various calculations, it is hard to take into account when using pure numerical simulation methodswang2018uniaxial .

Deformation of emergent cyrstals is essentially different from that of ordinary crystals: consider SkX as an example, its deformation is realized microscopically by variation of magnetization distribution, while a variation of the distance between two neighboring spins is not necessary. Based on this fundamental difference, emergent crystals possess two unique features compared with ordinary crystals: $\mathbf{a)}$ the local field configuration of the emergent particles inside an emergent crystal is geometrically unrestricted. The emergent particles in an emergent crystal are composed of multiple atoms or their electric and magnetic degrees of freedom with a nonconservative number. I.e., the emergent deformation is realized by absorbing the field configuration of neighboring atoms inside the emergent particle, for which its magnitude of deformability is unlimited. $\mathbf{b)}$ The emergent crystals can undergo emergent rotation with respect to the underlying material. We find that these two features dominates the nonelastic behavior of SkX suffering all kinds of external fields, including not only mechanical loading, but also magnetic fieldwang2017enhanced , electric fieldjiang2017direct , and temperature fieldmorikawa2017deformation . Moreover, these features permit a wide variety of possible crystalline structures of emergent crystals suffering different anisotropic fields, as shown in the results obtained here for SkX.

Strain-magnetic field Phase diagrams for MnSi

We further study the phase diagrams for MnSi suffering uniaxial distortion. We plot two types of $\varepsilon_{11}-b$ phase diagrams at the condition $t=0.5$, $\widetilde{A}_{e}=-0.05$: $\varepsilon_{11}-b$ phase diagram of equilibrium states (FIG. 4(a)), and $\varepsilon_{11}-b$ phase diagram of metastability for SkX (FIG. 4(b)). It is found that thermodynamically stable SkX exists in a very narrow range of $\varepsilon_{11}$, and exists in only one state: the RT-SkX phase. On the other hand, concerning the robustness of SkXkarube2016robust ; hu2018effect , metastable SkX exists in a much wider range of $\varepsilon_{11}$ and presence in a variety of different phases.

Concerning the wave nature of emergent crystalshu2018effect , all metastable states are stabilized by nonlinear mode-mode interaction, while the corresponding thermodynamically stable state: the generalized-conical (G-conical) phase is a single-Q structure whose wave vector is free to rotate in space. That is to say, the phase transition from any metastable emergent crystalline state to the G-conical is realized by cancelling of all the nonlinear mode-mode interactions (i.e., to set all related Fourier amplitudes to zero). This is physically difficult, because unlike annihilating a localized skyrmion, any Fourier component of magnetization exist globally and is hard to cancel through any local fluctuation. As a result, besides the well known topological protection attributed to isolated skyrmions, SkX is additionally protected globally by mode-mode interactions. In this sense, once we have a stable strain-free SkX phase, the appearance of various types of metastable SkX predicted in our calculation is anticipated, especially at low temperature range. This viewpoint is supported by the successful observation of the DS-SkX phase in uniaxially stretched FeGeshibata2015large .

Methods

$Emergent\,elasticity\,for\,2D\,hexagonal\,emergent\,\\ crystals\,in\,helimagnets.$

To formulate the theory of emergent elasticityhu2016emergent for emergent crystals in helimagnets, we switch from eq. (1) to a more generalized rescaled free energy functionalhu2017unified

[TABLE]

where

[TABLE]

denotes the elastic energy density of materials with cubic symmetry, $\gamma_{ij}=2\varepsilon_{ij}$ denote the engineering shear strains, and

[TABLE]

denotes the magnetoelastic energy density, where

[TABLE]

Eq. (1) can be derived from eq. (5) by setting all components of elastic strains to zero except $\varepsilon_{11}$. For deformable hexagonal 2-D emergent crystals, the rescaled magnetization $\mathbf{m}$ can be described mathematically as a Fourier serieshu2018wave :

[TABLE]

It is convenient to expand $\mathbf{m}_{\mathbf{q}_{ij}}$ as $\mathbf{m}_{\mathbf{q}_{ij}}=\tilde{c}_{ij1}\mathbf{P}_{ij1}+\tilde{c}_{ij2}\mathbf{P}_{ij2}+\tilde{c}_{ij3}\mathbf{P}_{ij3}$, where $\tilde{c}_{ij1}=\tilde{c}^{re}_{ij1}+{\rm i}\tilde{c}^{im}_{ij1}$, $\tilde{c}_{ij2}=\tilde{c}^{re}_{ij2}+{\rm i}\tilde{c}^{im}_{ij2}$, $\tilde{c}_{ij3}=\tilde{c}^{re}_{ij3}+{\rm i}\tilde{c}^{im}_{ij3}$ are complex variables to be determined, and $\mathbf{P}_{ij1}=\frac{1}{\sqrt{2}s_{i}q}[-iq_{ijy},iq_{ijx},s_{i}q]^{T}$, $\mathbf{P}_{ij2}=\frac{1}{s_{i}q}[q_{ijx},q_{ijy},0]^{T}$, $\mathbf{P}_{ij3}=\frac{1}{\sqrt{2}s_{i}q}[iq_{ijy},-iq_{ijx},s_{i}q]^{T}$ with $\mathbf{q}_{ij}=[q_{ijx},q_{ijy}]^{T}$, $|q_{ij}|=s_{i}q$. Assume that the undeformed wave vectors to be $q_{11}^{0}=q[0,1]^{T}$, $q_{12}^{0}=q[-\frac{\sqrt{3}}{2},-\frac{1}{2}]^{T}$ for a hexagonal lattice, the deformed wave vectors are related to the emergent elastic strains and emergent rotation angle byhu2016emergent

[TABLE]

where $s=1+\varepsilon^{e}_{11}+\varepsilon^{e}_{22}+\varepsilon^{e}_{11}\varepsilon^{e}_{22}-(\varepsilon^{e}_{12})^{2}+(\omega^{e})^{2}$. The averaged rescaled free energy density for emergent crystalline states can be obtained by substituting eqs. (9, 10) into eq. (5) and taking the volume average $\bar{w}=\frac{1}{V}\int\widetilde{w}dV$. At small disturbance around a metastable state and neglecting possible periodic elastic fields, the averaged rescaled free energy density for emergent crystalline states can be expanded in a quadratic form of all the independent variables as

[TABLE]

where $\bar{w}_{u}$ denotes the undisturbed averaged rescaled Helmholtz free energy density, and

[TABLE]

[TABLE]

[TABLE]

[TABLE]

One should notice that the form of eq. (14) has been simplified after canceling the degrees of freedom related to the in-plane rigid translationhu2016emergent . The linear constitutive equations are derived from eq. (11) as:

[TABLE]

where

[TABLE]

[TABLE]

[TABLE]

denote work conjugates of $\bm{\varepsilon}^{ea}$, $\bm{\varepsilon}$, and $\mathbf{m}^{a}$, respectively. When the material is subject to a disturbance of elastic strains, we have $d\mathbf{B}^{a}=\mathbf{0}$ and $d\bm{\sigma}^{ea}=\mathbf{0}$. In this case, after deduction we have from eq. (16)

[TABLE]

where

[TABLE]

Eq. (20) provides the linear relation between the emergent elastic strains, emergent rotations and the elastic strains, and $\bm{\lambda}$ is called the crossover strain ratio (CSR) matrix.

One should notice that the compliance matrices defined in eq. (15) provide necessary information to determine the metastability of a magnetic phase. For a state obtained from free energy minimization based on Fourier representation, the positive definiteness of the following matrixhu2018wave guarantees the intrinsic stability of the state:

[TABLE]

When the positive definiteness of $\bm{\Phi}$ is destroyed, the eigenvector corresponding to the smallest eigenvalue (softest mode) of $\bm{\Phi}$ indicates the direction of evolution for the system from the current state.

$Free\,energy\,minimization\,based\,on\,Monte-Carlo\\ Simulation.$

We perform Monte-Carlo Simulation based on annealing algorithm to minimize the free energy of the system. To process, we discretize the free energy density functional shown in eq. (1) with 2-D grids in Cartesian coordinates. Periodic boundary conditions are used in the calculation. For given elastic strains, the energy density functional can be written as

[TABLE]

where $m_{x}$, $m_{y}$, $m_{z}$ denote the components of the rescaled magnetization $\mathbf{m}$.

We use the central difference method to calculate the value of partial derivatives

[TABLE]

where the values of $\Delta x$ and $\Delta y$ determine the error of the Monte Carlo simulation.

In order to calculate the emergent strain of SkX, $\Delta x$ and $\Delta y$ are also set as variables. In this case, we calculated the averaged free energy density for different values of $\Delta x$ and $\Delta y$ to obtain the smallest one. And corresponding values of $\Delta x$ and $\Delta y$ determine the emergent elastic strains of SkX

[TABLE]

This method is only applicable to moderate emergent deformation of SkX. When the difference between equilibrium values of $\Delta x$ and $\Delta y$ is too large, the precision can no longer be guaranteed, and a new set of grids is needed.