Magnetostriction of helimagnets in the skyrmion crystal phase

Magnetic skyrmion, a topologically protected spin whirl, has attracted great attention for its potential applications in spintronics devices since its experimental discovery in 2009 [1]. Initially found in a small pocket-like region in the temperature–magnetic field phase diagram of bulk MnSi [1], skyrmions are later widely observed in other B20 bulk materials [2–4], thin films [5–9] and nanostructures [10–14] with enhanced stability. As a magnetic texture, skyrmion is not only affected by application of electromagnetic fields such as bias magnetic field, electric current [15], and laser [16], but also sensitive to non-electromagnetic fields such as temperature gradient [17] and mechanical stresses [3, 18]. While the interaction between skyrmions and the electromagnetic fields or the temperature gradient are well-understood both theoretically and experimentally, the effect of mechanical stress on the skyrmion is not well addressed, even though several previous experiments have proven that mechanical stress is an effective method for controlling the formation and stability of skyrmions [3, 18–24] and even their dynamic properties [25]. More specifically, the ultra-low emergent elastic stiffness [26] of the skyrmion crystal (SkX) observed in FeGe thin films [3] and in MnSi [18] makes it convenient to manipulate with mechanical stresses. Concerning the tensor properties of mechanical stresses and strains, sophisticated manipulation of the properties and stability of magnetic skyrmions may be realized by changing the type of mechanical loading, for which a comprehensive understanding of the interaction mechanism between magnetic skyrmions and mechanical loads is highly required.

The interaction between skyrmions and mechanical loads stem from the intrinsic magnetoelastic coupling in chiral magnets, which in turn affects the mechanical properties of the materials. Before the experimental discovery of SkX, the magnetostriction of a prototype helimagnet MnSi has been studied [27]. It is later found that appearance of the SkX phase leads to a jump of the magnetostriction [28] as well as the elastic constants [29, 30] as functions of the bias magnetic field. The first step toward a thorough understanding of the emergent elastic properties of the SkX phase is to explain these magnetoelastic-related phenomena of the SkX phase quantitatively within a unified theoretical framework. Recently, we have developed such a thermodynamic model concerning magnetoelastic effects for B20 compounds [31]. It explains quantitatively both the jump in the elastic constants in the conical-skyrmion (or skyrmion-conical) phase transition as the bias field increases [29, 30, 32] and the variation of the temperature–magnetic field phase diagram under uniaxial compression [18]. Nevertheless, magnetostriction in the SkX phase, which is the basic magnetoelastic property of any magnets, is not well understood either theoretically and experimentally.

In this work, we apply a combination of magnetostriction measurements and theoretical calculations to the real helimagnet MnSi, to demonstrate directly that the magnetostriction of B20 helimagnets are able to be well described by our developed theoretical model. Specifically, in the present work, we improve our previously-developed theory by employing more precise descriptions for different magnetic phases.

The rescaled Helmholtz free energy density for the chiral magnets [31] in this study can be written as

$$\begin{eqnarray}&&\widetilde{w}\left({\bf{m}},{\varepsilon }_{{ij}}\right)={\widetilde{w}}_{m}({\bf{m}})+{\widetilde{w}}_{e}\left({\bf{m}},{\varepsilon }_{{ij}}\right),\end{eqnarray} \tag{ 1 }$$

where m denotes the rescaled magnetization, ε_ij denotes the elastic strains, and ${\widetilde{w}}_{m}({\bf{m}})$ denotes the part of the free energy density that depends solely on the magnetization. It can be written as

$$\begin{eqnarray}\begin{array}{rcl}{\widetilde{w}}_{m}({\bf{m}}) & = & \displaystyle \sum _{i=1}^{3}{\left(\displaystyle \frac{\partial {\bf{m}}}{\partial {r}_{i}}\right)}^{2}+2{\bf{m}}\cdot ({\boldsymbol{\nabla }}\times {\bf{m}})-2{\bf{b}}\cdot {\bf{m}}\\ & & +\ t{{\bf{m}}}^{2}+{{\bf{m}}}^{4}+\displaystyle \sum _{i=1}^{3}{\widetilde{A}}_{e}{\left(\displaystyle \frac{\partial {m}_{i}}{\partial {r}_{i}}\right)}^{2},\end{array}\end{eqnarray} \tag{ 2 }$$

Here, the terms on the r.h.s. of equation (2) denote, respectively, the exchange energy density, the Dzyaloshinskii-Moriya (DM) interaction, the Zeeman energy density with the rescaled magnetic field b, the second and fourth order Landau expansion terms with the rescaled temperature t, and the exchange anisotropy with a rescaled coefficient ${\widetilde{A}}_{e}$.

The second term on the r.h.s.of equation (1) denotes the part of the free energy density that depends on the elastic strains ${\varepsilon }_{{ij}};$ it can be expanded as

$$\begin{eqnarray}&&{\widetilde{w}}_{e}\left({\bf{m}},{\varepsilon }_{{ij}}\right)={\widetilde{w}}_{{el}}\left({\varepsilon }_{{ij}}\right)+{\widetilde{w}}_{{me}}\left({\bf{m}},{\varepsilon }_{{ij}}\right),\end{eqnarray} \tag{ 3 }$$

where ${\widetilde{w}}_{{el}}\left({\varepsilon }_{{ij}}\right)$ denotes the rescaled elastic energy density, and ${\widetilde{w}}_{{me}}\left({\bf{m}},{\varepsilon }_{{ij}}\right)$ denotes the rescaled magnetoelastic free energy density. Expressions for ${\widetilde{w}}_{{el}}\left({\varepsilon }_{{ij}}\right)$ and ${\widetilde{w}}_{{me}}\left({\bf{m}},{\varepsilon }_{{ij}}\right)$, as well as the rescaling process of equation (2), are given in the appendix.

2.1. Magnetostriction strains and the temperature–magnetic field phase diagram calculated for different magnetic phases suffering mechanical loads

Magnetostriction strains refer to the homogeneous elastic strains caused by magnetization when no mechanical load is applied. As a result, they can be derived from $\tfrac{\partial \overline{w}\left({\bf{m}},{\varepsilon }_{{ij}}\right)}{\partial {\varepsilon }_{{ij}}}=0$ when i = j, and $\tfrac{\partial \overline{w}\left({\bf{m}},{\varepsilon }_{{ij}}\right)}{\partial {\gamma }_{{ij}}}=0$ when $i\ne j$, where $\overline{w}\left({\bf{m}},{\varepsilon }_{{ij}}\right)=\tfrac{1}{V}\int \widetilde{w}\left({\bf{m}},{\varepsilon }_{{ij}}\right){\rm{d}}{V}$. The analytical solution of the magnetostriction strains can be obtained from ${\varepsilon }_{{ij}}=\tfrac{1}{V}\int {\varepsilon }_{{ij}}^{* }{\rm{d}}{V}$ with

$$\begin{eqnarray}\begin{array}{rcl}{\varepsilon }_{11}^{* } & = & {K}^{* }{m}^{2}-{L}_{1}^{* }{m}_{1}^{2}-{L}_{2}^{* }{m}_{3}^{2}+{L}_{O1}^{* }\left({m}_{3}{m}_{\mathrm{1,2}}-{m}_{2}{m}_{\mathrm{1,3}}\right)+{L}_{O2}^{* }\left({m}_{3}{m}_{\mathrm{2,1}}-{m}_{2}{m}_{\mathrm{3,1}}\right)\\ & & +\ {L}_{O3}^{* }{m}_{1}\left({m}_{\mathrm{2,3}}-{m}_{\mathrm{3,2}}\right),\\ {\varepsilon }_{22}^{* } & = & {K}^{* }{m}^{2}-{L}_{1}^{* }{m}_{2}^{2}-{L}_{2}^{* }{m}_{1}^{2}+{L}_{O1}^{* }\left({m}_{1}{m}_{\mathrm{2,3}}-{m}_{3}{m}_{\mathrm{2,1}}\right)+{L}_{O2}^{* }\left({m}_{1}{m}_{\mathrm{3,2}}-{m}_{3}{m}_{\mathrm{1,2}}\right)\\ & & +\ {L}_{O3}^{* }{m}_{2}\left({m}_{\mathrm{3,1}}-{m}_{\mathrm{1,3}}\right),\\ {\varepsilon }_{33}^{* } & = & {K}^{* }{m}^{2}-{L}_{1}^{* }{m}_{3}^{2}-{L}_{2}^{* }{m}_{2}^{2}+{L}_{O1}^{* }\left({m}_{2}{m}_{\mathrm{3,1}}-{m}_{1}{m}_{\mathrm{3,2}}\right)+{L}_{O2}^{* }\left({m}_{2}{m}_{\mathrm{1,3}}-{m}_{1}{m}_{\mathrm{2,3}}\right)\\ & & +\ {L}_{O3}^{* }{M}_{3}\left({m}_{\mathrm{1,2}}-{m}_{\mathrm{2,1}}\right),\\ {\gamma }_{23}^{* } & = & {C}_{44}^{-1}[-{\tilde{L}}_{3}{m}_{2}{m}_{3}+{\tilde{L}}_{O6}{m}_{1}({m}_{\mathrm{2,2}}-{m}_{\mathrm{3,3}})+{m}_{2}({\tilde{L}}_{O4}{m}_{\mathrm{1,2}}+{\tilde{L}}_{O5}{m}_{\mathrm{2,1}})\\ & & -\ {m}_{3}({\tilde{L}}_{O4}{m}_{\mathrm{1,3}}+{\tilde{L}}_{O5}{m}_{\mathrm{3,1}})],\\ {\gamma }_{13}^{* } & = & {C}_{44}^{-1}[-{\tilde{L}}_{3}{m}_{1}{m}_{3}+{\tilde{L}}_{O6}{m}_{2}({m}_{\mathrm{3,3}}-{m}_{\mathrm{1,1}})+{m}_{3}({\tilde{L}}_{O4}{m}_{\mathrm{2,3}}+{\tilde{L}}_{O5}{m}_{\mathrm{3,2}})\\ & & -\ {m}_{1}({\tilde{L}}_{O4}{m}_{\mathrm{2,1}}+{\tilde{L}}_{O5}{m}_{\mathrm{1,2}})],\\ {\gamma }_{12}^{* } & = & {C}_{44}^{-1}[-{\tilde{L}}_{3}{m}_{1}{m}_{2}+{\tilde{L}}_{O6}{m}_{3}({m}_{\mathrm{1,1}}-{m}_{\mathrm{2,2}})+{m}_{1}({\tilde{L}}_{O4}{m}_{\mathrm{3,1}}+{\tilde{L}}_{O5}{m}_{\mathrm{1,3}})\\ & & -\ {m}_{2}({\tilde{L}}_{O4}{m}_{\mathrm{3,2}}+{\tilde{L}}_{O5}{m}_{\mathrm{2,3}})],\end{array}\end{eqnarray} \tag{ 4 }$$

In equation (4), ${\gamma }_{{ij}}=2{\varepsilon }_{{ij}},i\ne j$ denotes the engineering shear strains. The rescaled magnetization ${\bf{m}}={\left[{m}_{1},{m}_{2},{m}_{3}\right]}^{T}$, ${m}^{2}={m}_{1}^{2}+{m}_{2}^{2}+{m}_{3}^{2}$ and m_i,j (i, j = 1, 2, 3) is the first-order partial derivatives. K*, ${L}_{i}^{* }$ (i = 1, 2, 3), and ${L}_{{Oi}}^{* }$ (i = 1, 2, ⋯, 6) are coefficients of different orders of magnetoelastic interactions and are explained in detail in the appendix. Substitution of the expression of rescaled magnetization m for different magnetic phases into equation (4) yields detailed expression for the magnetostriction strains. To determine the exact values of magnetostriction strains at a given temperature and magnetic field, one first minimizes $\overline{w}\left({\bf{m}},{\varepsilon }_{{ij}}({\bf{m}}\right)$ with respect to m for different magnetic phases at the condition considered [31]. By changing the values of temperature and magnetic field and then repeating this process, we obtain the temperature–magnetic field phase diagram.

We choose the prototype chiral magnet MnSi as an example to calculate the temperature–magnetic field phase diagram using the thermodynamic parameters listed in⁴ . To stabilize the SkX phase, we introduce in the free energy functional a negative exchange anisotropy [33], and describe the SkX phase within the fourth order Fourier representation [34] instead of the third-order Fourier representation used in our previous work [31], from which we obtain a lower free energy density for the SkX phase. For the conical and helical phases, we include a 'double-Q' Fourier term in the description of the magnetization with the same direction of the wave vector [35]. From this, we obtain a lower free energy densities for the conical and helical phases. With these improvements, both the calculated temperature–magnetic field phase diagram and the magnetostriction strain agree better with the corresponding experiments.

The calculated results reproduce the typical features of chiral magnets, as shown in figure 1. Below the critical temperature (T_c ∼ 30 K for MnSi),we obtain a helimagnetic ground state due to the competition between the exchange and DM interactions. The helimagnetic state transforms into a conical phase and finally into the ferromagnetic state with increasing magnetic field. Skyrmions, in the form of the crystal phase, appear very close to T_c only in a tiny region of magnetic field and temperature.

Zoom In Zoom Out Reset image size

Using the values of the equilibrium magnetization obtained from these phase diagram calculations, we can further calculate the magnetostriction strains at any given temperature and magnetic field through equation (4). Here, we chose the variation of Δε₃₃ as an example, because it is easily measured experimentally. This quantity, Δε₃₃ = ε₃₃(b) − ε₃₃(0), denotes the change of ε₃₃ due to a change of the magnetic field.

The typical magnetic field dependence of Δε₃₃ at the temperature T = 26.4 K is shown in figure 2(a). As the field increases from zero, the magnetostriction first decreases up to the lower critical field B_C1 ≈ 0.07 T. It then increases approximately quadratically with the field strengthen up to the upper critical magnetic field B_C2 ≈ 0.25 T. This is accompanied by a hill-like jump between the two fields B_A1 ≈ 0.08 T and B_A2 ≈ 0.11 T. At magnetic field higher than B > B_C2, we obtain an approximately linear dependence of the magnetostriction on the field strength. Compared with the conventional magnetic phases in the prototypical helimagnet MnSi and in previous work [27], B_C1 and B_C2 are easily identified to correspond the critical fields for the helical-to-conical phase and the conical-to-ferromagnetic phase transition, respectively. The hill-like jump in the field range [B_A1, B_A2] characterizes the skyrmion phase.

Zoom In Zoom Out Reset image size

The detailed dependence of the magnetostriction on the magnetic field at various temperatures below T_c is shown in figure 2(b). This jumped behavior emerges in the whole region of the skyrmion phase.

3.1. Derivation of the analytical expression of magnetostriction strains in bulk MnSi

To explain this jump, we derive below analytical expression of magnetostriction strains for the SkX phase within the triple-Q approximation of the magnetization

$$\begin{eqnarray}&&\begin{array}{c}\begin{array}{c}\begin{array}{c}\begin{array}{c}\begin{array}{c}\begin{array}{c}\begin{array}{c}{{\bf{m}}}_{{\rm{t}}{\rm{r}}{\rm{i}}{\rm{p}}{\rm{l}}{\rm{e}}Q}=\left[\begin{array}{c}0\\ 0\\ {m}_{0}\end{array}\right]+\displaystyle \frac{\sqrt{3}{m}_{q}}{3}\left\{\left[\begin{array}{c}0\\ \sin ({{\bf{q}}}_{1}{\bf{r}})\\ -\cos ({{\bf{q}}}_{1}{\bf{r}})\end{array}\Space{0ex}{2.0ex}{0ex}]\right.\right.+\left[\begin{array}{c}-\displaystyle \frac{\sqrt{3}}{2}\sin ({{\bf{q}}}_{2}{\bf{r}})\\ -\displaystyle \frac{1}{2}\sin ({{\bf{q}}}_{2}{\bf{r}})\\ -\cos ({{\bf{q}}}_{2}{\bf{r}})\end{array}\right]+\left.\Space{0ex}{2.0ex}{0ex}\begin{array}{c}\displaystyle \frac{\sqrt{3}}{2}\sin ({{\bf{q}}}_{3}{\bf{r}})\\ -\displaystyle \frac{1}{2}\sin ({{\bf{q}}}_{3}{\bf{r}})\\ -\cos ({{\bf{q}}}_{3}{\bf{r}})\end{array}\Space{0ex}{2.0ex}{0ex}]\right\},\end{array}\end{array}\end{array}\end{array}\end{array}\end{array}\end{array}\end{eqnarray} \tag{ 5 }$$

with ${{\bf{q}}}_{1}=q{\left[\begin{array}{lll}1 & 0 & 0\end{array}\right]}^{T}$, ${{\bf{q}}}_{2}=q{\left[\begin{array}{lll}-\tfrac{1}{2} & \tfrac{\sqrt{3}}{2} & 0\end{array}\right]}^{{T}}$, ${{\bf{q}}}_{3}=q{\left[\begin{array}{lll}-\tfrac{1}{2} & -\tfrac{\sqrt{3}}{2} & 0\end{array}\right]}^{T}$, which yields from equation (4)

$$\begin{eqnarray}\begin{array}{rcl}{\left({\varepsilon }_{11}\right)}_{{\rm{skyrmion}}} & = & \left[{K}^{* }\left({m}_{q}^{2}+{m}_{0}^{2}\right)+\displaystyle \frac{1}{4}\left(-{L}_{1}^{* }-2{L}_{2}^{* }+{L}_{O1}^{* }q-2{L}_{O2}^{* }q+{L}_{O3}^{* }q\right){m}_{q}^{2}-{L}_{2}^{* }{m}_{0}^{2}\right],\\ {\left({\varepsilon }_{22}\right)}_{{\rm{skyrmion}}} & = & \left[{K}^{* }\left({m}_{q}^{2}+{m}_{0}^{2}\right)+\displaystyle \frac{1}{4}\left(-{L}_{1}^{* }-{L}_{2}^{* }+{L}_{O1}^{* }q-2{L}_{O2}^{* }q+{L}_{O3}^{* }q\right){m}_{q}^{2}\right],\\ {\left({\varepsilon }_{33}\right)}_{{\rm{skyrmion}}} & = & \left[{K}^{* }\left({m}_{q}^{2}+{m}_{0}^{2}\right)+\displaystyle \frac{1}{4}\left(-2{L}_{1}^{* }-{L}_{2}^{* }+2{L}_{O1}^{* }q+2{L}_{O3}^{* }q\right){m}_{q}^{2}-{L}_{1}^{* }{m}_{0}^{2}\right],\\ {\bar{\gamma }}_{{ij}} & = & 0\ {\rm{when}}\ i\ne j.\end{array}\end{eqnarray} \tag{ 6 }$$

For the conical phase, we have ${{\bf{m}}}_{{\rm{conical}}}={\left[\begin{array}{lll}{m}_{{qc}}\cos \left({{qr}}_{3}\right) & {m}_{{qc}}\sin \left({{qr}}_{3}\right) & {m}_{0c}\end{array}\right]}^{T}$,which gives from equation (4)

$$\begin{eqnarray}\begin{array}{c}\begin{array}{rcl}{\left({\varepsilon }_{11}\right)}_{{\rm{c}}{\rm{o}}{\rm{n}}{\rm{i}}{\rm{c}}{\rm{a}}{\rm{l}}} & = & \left[{K}^{\ast }\left({m}_{qc}^{2}+{m}_{0c}^{2}\right)+\displaystyle \frac{1}{2}(-{L}_{1}^{\ast }+{L}_{O1}^{\ast }q+{L}_{O3}^{\ast }q){m}_{qc}^{2}-{L}_{2}^{\ast }{m}_{0c}^{2}\right],\\ {\left({\varepsilon }_{22}\right)}_{{\rm{c}}{\rm{o}}{\rm{n}}{\rm{i}}{\rm{c}}{\rm{a}}{\rm{l}}} & = & \left[{K}^{\ast }\left({m}_{qc}^{2}+{m}_{0c}^{2}\right)+\displaystyle \frac{1}{2}(-{L}_{1}^{\ast }+{L}_{O1}^{\ast }q+{L}_{O3}^{\ast }q){m}_{qc}^{2}-\displaystyle \frac{1}{2}{L}_{2}^{\ast }{m}_{qc}^{2}\right],\\ {\left({\varepsilon }_{33}\right)}_{{\rm{c}}{\rm{o}}{\rm{n}}{\rm{i}}{\rm{c}}{\rm{a}}{\rm{l}}} & = & \left[{K}^{\ast }\left({m}_{qc}^{2}-{m}_{0c}^{2}\right)+{L}_{1}^{\ast }{m}_{0c}^{2}-{L}_{O2}^{\ast }{qm}_{qc}^{2}-\displaystyle \frac{1}{2}{L}_{2}^{\ast }{m}_{qc}^{2}\right],\\ {\bar{\gamma }}_{ij} & = & 0\,{\rm{w}}{\rm{h}}{\rm{e}}{\rm{n}}\,i\ne j.\end{array}\end{array}\end{eqnarray} \tag{ 7 }$$

Comparing equations (6) and (7), we have

$$\begin{eqnarray}&&\begin{array}{l}{\left({\varepsilon }_{11}\right)}_{{\rm{skyrmion}}}-{\left({\varepsilon }_{11}\right)}_{{\rm{conical}}}\approx \displaystyle \frac{1}{4}({L}_{1}^{* }-2{L}_{2}^{* }){m}_{q}^{2},\\ {\left({\varepsilon }_{22}\right)}_{{\rm{skyrmion}}}-{\left({\varepsilon }_{22}\right)}_{{\rm{conical}}}\approx \displaystyle \frac{1}{4}({L}_{2}^{* }+{L}_{1}^{* }){m}_{q}^{2},\\ {\left({\varepsilon }_{33}\right)}_{{\rm{skyrmion}}}-{\left({\varepsilon }_{33}\right)}_{{\rm{conical}}}\approx \displaystyle \frac{1}{4}({L}_{2}^{* }-2{L}_{1}^{* }){m}_{q}^{2},\end{array}\end{eqnarray} \tag{ 8 }$$

where it is found through calculation that m₀ ≈ m_0c and m_q ≈ m_qc, and the condition $\left|{L}_{O1}^{* }q\right|$, $\left|{L}_{O2}^{* }q\right|$, $\left|{L}_{O3}^{* }q\right|\ll \left|{L}_{1}^{* }\right|$, $\left|{L}_{2}^{* }\right|$ is used. The magnitude of the jump in ε₃₃ through the phase transition from the conical phase to the SkX phase is determined by the coefficient ${L}_{2}^{* }-2{L}_{1}^{* }$. For MnSi ${L}_{2}^{* }-2{L}_{1}^{* }\gt 0$, which yields the result in figure 2. Another phenomenon shown in figure 2(b) is that the jump in ε₃₃ in the conical-skyrmion phase transition becomes more insignificant as the temperature rises. As shown in equation (6), the difference in ε₃₃ between the SkX phase and the conical phase depends not only on the parameter ${L}_{2}^{* }-2{L}_{1}^{* }$ but also on the value of m_q. Since the skyrmion phase appears near the Curie temperature of the material, and near the Curie temperature we have ${m}_{q}\to 0$, the difference of ε₃₃ in the two phases gradually vanishes as the temperature rises. According to equation (7), this phenomenon is generally true for all strain components.

To test the theoretical results, we measured the magnetostriction of MnSi using the capacitance method [36]. We used a rectangular, single-crystal sample of MnSi that was 1.82 mm long, 1.09 mm wide, and 2.61 mm high, with the pair of end-faces parallel to the (001) crystallographic plane. We always measured the sample length in the [100] direction, which is also parallel to the orientation of the magnetic field.

The complete variation of Δε₃₃ (calculated from ΔL/L, where L denotes the length in the height direction and ΔL denotes its variation) with the applied magnetic field for MnSi at different temperatures is shown in figure 3(a). The rules for defining the four critical points in figure 3(a) are as follows: B_C1 denotes the field strength at which the first-order derivative dB/dT is zero; this corresponds to the minimum field in the low magnetic field region, which is marked by the black arrow in figure 3(b). The critical points B_A1 and B_A2 are defined as the maximum and minimum of the first-order derivative curve in the intermediate magnetic field region, as shown in figure 3(b). These two points define the region where skyrmions exist. The last transition field corresponds to an abrupt change in slope, as obviously shown by the sudden jump of ΔL/L in the derivative curve. This critical field B_C2 in figure 3(a) is determined from the midpoint of the transition in d(ΔL/L)/dB, as marked by the black arrow in the derivative curve shown in figure 3(b). Figure 3(c) shows the experimental phase diagram for MnSi, which is constructed from figure 3(a).

Zoom In Zoom Out Reset image size

We have shown that the experimental measurements well verified the theoretical results. The magnitude of the maximum jumped value of Δε₃₃ is on the order of 2.5 × 10⁻⁷, the same as the theoretical result. Nevertheless, three discrepancies exist between experiment and theory. The first is that the jump in Δε₃₃ between the magnetic phases is smoother in the experiment. The reason for this discrepancy is probably due to neglect of possible mixed states of the conical and SkX phases in the theoretical calculations. As shown by various TEM(Transmission Electron Microscopy) observations of the SkX phase [5, 6], at the critical magnetic field for the conical-SkX transition, some isolated skyrmions or skyrmion clusters appear due to the presence of defects or grain boundaries, thus forming a mixed state of conical and skyrmions phases, which smooths the jump in Δε₃₃. The second concern is the intermediate magnetic (IM) state, which appears in a small window between the paramagnetic phase and the skyrmion phase in figure 3(c), but which is not present in the theoretical phase diagram shown in figure 1. In the IM phase, the magnetostriction strain changes smoothly, as illustrated by the curve for Δε₃₃ at T = 29.6 K in figure 3(a), which indicates that the helical-conical transition is missing. This IM state may be related to the precursor states discussed in previous works [37–41] , but further exploration is needed to confirm this. The third concern is the phase transition at high temperature from calculated phase diagram is not consistent with the experimental one. The SkX phase is connecting to the helical phase at high temperature in theoretical calculation, but in experiment the SX phase is connecting to conical phase. This may be due to the fact that theoretical calculation is performed within mean-field approximation, where the effect of fluctuation is neglected.

In summary, we have derived a general expression of magnetostriction strains for B20 helimagnets in the SkX phase. For MnSi, we calculated the normal strain ε₃₃ at different temperatures and magnetic fields, and we also perform the corresponding experimental measurements, which agrees quantitatively with each other. ε₃₃ undergoes a sudden jump that accompany a conical-skyrmion phase transition at any temperature, which can be used to confirm the presence of the SkX phase.

This work was supported by the National Key R&D Program of China, grant no.2017YFA0303201, the Natural Science Foundation of China, grant no. 11772360, 51622105, the Key Research Program of Frontier Sciences, CAS, grant no. QYZDB-SSWSLH009, the Key Research Program of the Chinese Academy of Sciences, grant no. KJZDSWM01,Pearl River Nova Program of Guangzhou, grant no.201806010134.

The free energy density functional for cubic helimagnets incorporating magnetoelastic interactions can be written as

$$\begin{eqnarray}\begin{array}{rcl}w\left({\bf{M}},{\varepsilon }_{{ij}}\right) & = & \displaystyle \sum _{i=1}^{3}A{\left(\displaystyle \frac{\partial {\bf{M}}}{\partial {x}_{i}}\right)}^{2}+D{\bf{M}}\cdot ({\boldsymbol{\nabla }}\times {\bf{M}})-{\bf{B}}\cdot {\bf{M}}+\alpha \left(T-{T}_{0}\right){{\bf{M}}}^{2}+\beta {{\bf{M}}}^{4}\\ & & +\ \displaystyle \sum _{i=1}^{3}{A}_{e}{\left(\displaystyle \frac{\partial {M}_{i}}{\partial {x}_{i}}\right)}^{2}+{w}_{{el}}+{w}_{{me}},\end{array}\end{eqnarray} \tag{ A.1 }$$

where

$$\begin{eqnarray}&&{w}_{{el}}=\displaystyle \frac{1}{2}{C}_{11}\left({\varepsilon }_{11}^{2}+{\varepsilon }_{22}^{2}+{\varepsilon }_{33}^{2}\right)+{C}_{12}({\varepsilon }_{11}{\varepsilon }_{22}+{\varepsilon }_{11}{\varepsilon }_{33}+{\varepsilon }_{22}{\varepsilon }_{33})+\displaystyle \frac{1}{2}{C}_{44}\left({\gamma }_{12}^{2}+{\gamma }_{13}^{2}+{\gamma }_{23}^{2}\right),\end{eqnarray} \tag{ A.2 }$$

$$\begin{eqnarray}\begin{array}{rcl}{w}_{{me}} & = & \displaystyle \frac{1}{{M}_{S}^{2}}[{{KM}}^{2}{\varepsilon }_{{ii}}+{L}_{1}\left({M}_{1}^{2}{\varepsilon }_{11}+{M}_{2}^{2}{\varepsilon }_{22}+{M}_{3}^{2}{\varepsilon }_{33}\right)+{L}_{2}\left({M}_{3}^{2}{\varepsilon }_{11}+{M}_{1}^{2}{\varepsilon }_{22}+{M}_{2}^{2}{\varepsilon }_{33}\right)\\ & & +\ {L}_{3}({M}_{1}{M}_{2}{\gamma }_{12}+{M}_{1}{M}_{3}{\gamma }_{13}+{M}_{2}{M}_{3}{\gamma }_{23})]+\displaystyle \frac{1}{{M}_{s}^{2}}\displaystyle \sum _{i=1}^{6}{L}_{{Oi}}{f}_{{Oi}},\end{array}\end{eqnarray} \tag{ A.3 }$$

and

$$\begin{eqnarray}\begin{array}{rcl}{f}_{O1} & = & {\varepsilon }_{11}\left({M}_{\mathrm{1,2}}{M}_{3}-{M}_{\mathrm{1,3}}{M}_{2}\right)+{\varepsilon }_{22}\left({M}_{\mathrm{2,3}}{M}_{1}-{M}_{\mathrm{2,1}}{M}_{3}\right)+{\varepsilon }_{33}\left({M}_{\mathrm{3,1}}{M}_{2}-{M}_{\mathrm{3,2}}{M}_{1}\right),\\ {f}_{O2} & = & {\varepsilon }_{11}\left({M}_{\mathrm{3,1}}{M}_{2}-{M}_{\mathrm{2,1}}{M}_{3}\right)+{\varepsilon }_{22}\left({M}_{\mathrm{1,2}}{M}_{3}-{M}_{\mathrm{3,2}}{M}_{1}\right)+{\varepsilon }_{33}\left({M}_{\mathrm{2,3}}{M}_{1}-{M}_{\mathrm{1,3}}{M}_{2}\right),\\ {f}_{O3} & = & {\varepsilon }_{11}{M}_{1}\left({M}_{\mathrm{2,3}}-{M}_{\mathrm{3,2}}\right)+{\varepsilon }_{22}{M}_{2}\left({M}_{\mathrm{3,1}}-{M}_{\mathrm{1,3}}\right)+{\varepsilon }_{33}{M}_{3}\left({M}_{\mathrm{1,2}}-{M}_{\mathrm{2,1}}\right),\\ {f}_{O4} & = & {\gamma }_{23}\left({M}_{\mathrm{1,3}}{M}_{3}-{M}_{\mathrm{1,2}}{M}_{2}\right)+{\gamma }_{13}\left({M}_{\mathrm{2,1}}{M}_{1}-{M}_{\mathrm{2,3}}{M}_{3}\right)+{\gamma }_{12}\left({M}_{\mathrm{3,2}}{M}_{2}-{M}_{\mathrm{3,1}}{M}_{1}\right),\\ {f}_{O5} & = & {\gamma }_{23}\left({M}_{\mathrm{3,1}}{M}_{3}-{M}_{\mathrm{2,1}}{M}_{2}\right)+{\gamma }_{13}\left({M}_{\mathrm{1,2}}{M}_{1}-{M}_{\mathrm{3,2}}{M}_{3}\right)+{\gamma }_{12}\left({M}_{\mathrm{2,3}}{M}_{2}-{M}_{\mathrm{1,3}}{M}_{1}\right),\\ {f}_{O6} & = & {\gamma }_{23}{M}_{1}\left({M}_{\mathrm{3,3}}-{M}_{\mathrm{2,2}}\right)+{\gamma }_{13}{M}_{2}\left({M}_{\mathrm{1,1}}-{M}_{\mathrm{3,3}}\right)+{\gamma }_{12}{M}_{3}\left({M}_{\mathrm{2,2}}-{M}_{\mathrm{1,1}}\right).\end{array}\end{eqnarray} \tag{ A.4 }$$

Here ${\bf{M}}={\left[{M}_{1},{M}_{2},{M}_{3}\right]}^{T}$ denotes the magnetization vector, M_s denotes the saturation magnetization, and ${M}^{2}={M}_{1}^{2}+{M}_{2}^{2}+{M}_{3}^{2}$. The first term in equation (A.1) describes the exchange energy density with stiffness A; the second term is the Dzyaloshinskii–Moriya (DM) interaction with coupling constant D, which determines the period and direction of the periodic magnetization; the third term is the Zeeman energy density in the applied magnetic field B; $\alpha \left(T-{T}_{0}\right){{\bf{M}}}^{2}+\beta {{\bf{M}}}^{4}$ are two Landau expansion terms; A_e denotes the coefficient of the exchange anisotropy; and w_el and w_me denote respectively the elastic energy density given in equation (A.2) and magnetoelastic free energy density given in equation (A.3).

Equation (A.1) provides an implicit model for studying the effects of magnetoelastic interactions on the SkX phase. The effects of these interactions cannot be understood by simply examining the magnetoelastic thermodynamic parameters, such as K, L₁, etc., because they are also related to magnetic thermodynamic parameters, such as A and D. It is thus more convenient to write the free energy density functional in a rescaled form given in equation (A.1), where

$$\begin{eqnarray}&&\widetilde{w}\left({\bf{m}},{\varepsilon }_{{ij}}\right)=\displaystyle \frac{\beta }{{G}^{2}}w\left({\bf{M}},{\varepsilon }_{{ij}}\right).\end{eqnarray} \tag{ A.5 }$$

$$\begin{eqnarray}&&{\widetilde{w}}_{{el}}=\displaystyle \frac{1}{2}{\tilde{C}}_{11}\left({\varepsilon }_{11}^{2}+{\varepsilon }_{22}^{2}+{\varepsilon }_{33}^{2}\right)+{\tilde{C}}_{12}({\varepsilon }_{11}{\varepsilon }_{22}+{\varepsilon }_{11}{\varepsilon }_{33}+{\varepsilon }_{22}{\varepsilon }_{33})+\displaystyle \frac{1}{2}{\tilde{C}}_{44}\left({\gamma }_{12}^{2}+{\gamma }_{13}^{2}+{\gamma }_{23}^{2}\right).\end{eqnarray} \tag{ A.6 }$$

$$\begin{eqnarray}&&\begin{array}{r}{\tilde{w}}_{{me}}=\tilde{K}{m}^{2}{\varepsilon }_{{ii}}+{\tilde{L}}_{1}({m}_{1}^{2}{\varepsilon }_{11}+{m}_{2}^{2}{\varepsilon }_{22}+{m}_{3}^{2}{\varepsilon }_{33})+{\tilde{L}}_{2}({m}_{3}^{2}{\varepsilon }_{11}+{m}_{1}^{2}{\varepsilon }_{22}+{m}_{2}^{2}{\varepsilon }_{33})\\ \,\,\,+\,{\tilde{L}}_{3}({m}_{1}{m}_{2}{\gamma }_{12}+{m}_{1}{m}_{3}{\gamma }_{13}+{m}_{2}{m}_{3}{\gamma }_{23})+\displaystyle \sum _{i=1}^{6}{\tilde{L}}_{{Oi}}{\tilde{f}}_{{Oi}}.\end{array}\end{eqnarray} \tag{ A.7 }$$

and

$$\begin{eqnarray}\begin{array}{rcl}{\bf{r}} & = & \displaystyle \frac{{\bf{x}}}{{L}_{D}},{\bf{b}}=\displaystyle \frac{{\bf{B}}}{B},{\bf{m}}=\displaystyle \frac{{\bf{M}}}{{M}_{0}},{L}_{D}=\displaystyle \frac{2A}{D},G=\displaystyle \frac{{D}^{2}}{4A},\\ B & = & 2{{GM}}_{0},{M}_{0}=\sqrt{\displaystyle \frac{G}{\beta }},t=\displaystyle \frac{\alpha \left(T-{T}_{0}\right)}{G}.\end{array}\end{eqnarray} \tag{ A.8 }$$

The rescaled thermodynamic parameters (parameters with wavy overlines) are defined by

$$\begin{eqnarray}\begin{array}{rcl}{\tilde{A}}_{e} & = & \displaystyle \frac{{A}_{e}}{A},\,\tilde{K}=\displaystyle \frac{K}{{GM}_{s}^{2}},\,{\tilde{L}}_{1}=\displaystyle \frac{{L}_{1}}{{GM}_{S}^{2}},\,{\tilde{L}}_{2}=\displaystyle \frac{{L}_{2}}{{GM}_{s}^{2}},\,{\tilde{L}}_{3}=\displaystyle \frac{{L}_{3}}{{GM}_{s}^{2}},\\ & = & {\tilde{L}}_{Oi}=\displaystyle \frac{2{L}_{Oi}}{{DM}_{S}^{2}},(i=1,2,\,\ldots ,\,6),\,{\tilde{C}}_{ij}=\displaystyle \frac{\beta }{{G}^{2}}{C}_{ij}.\end{array}\end{eqnarray} \tag{ A.9 }$$

In equation (A.7), ${\tilde{f}}_{{Oi}}$, (i = 1, 2, ..., 6) can be derived from equation (A.4) by replacing M by m. In equation (7) of the main text, coefficients with the superscript '*' are defined by

$$\begin{eqnarray*}\begin{array}{rcl}{K}^{* } & = & \displaystyle \frac{-{\tilde{C}}_{11}\widetilde{K}+{\tilde{C}}_{12}\left(\widetilde{K}+{\widetilde{L}}_{1}+{\tilde{L}}_{2}\right)}{\left({\tilde{C}}_{11}-{\tilde{C}}_{12}\right)\left({\tilde{C}}_{11}+2{\tilde{C}}_{12}\right)},\\ {L}_{1}^{* } & = & \displaystyle \frac{{\tilde{L}}_{1}}{\left({\tilde{C}}_{11}-{\tilde{C}}_{12}\right)},{L}_{2}^{* }=\displaystyle \frac{{\tilde{L}}_{2}}{\left({\tilde{C}}_{11}-{\tilde{C}}_{12}\right)},\\ {L}_{O1}^{* } & = & \displaystyle \frac{-{\tilde{C}}_{11}{\tilde{L}}_{O1}+{\tilde{C}}_{12}\left(-{\tilde{L}}_{O1}+{\tilde{L}}_{O2}+{\tilde{L}}_{O3}\right)}{\left({\tilde{C}}_{11}-{\tilde{C}}_{12}\right)\left({\tilde{C}}_{11}+2{\tilde{C}}_{12}\right)},\\ {L}_{O2}^{* } & = & \displaystyle \frac{{\tilde{C}}_{11}{\tilde{L}}_{O2}-{\tilde{C}}_{12}\left({\tilde{L}}_{O1}-{\tilde{L}}_{O2}+{\tilde{L}}_{O3}\right)}{\left({\tilde{C}}_{11}-{\tilde{C}}_{12}\right)\left({\tilde{C}}_{11}+2{\tilde{C}}_{12}\right)},\\ {L}_{O3}^{* } & = & \displaystyle \frac{{\tilde{C}}_{12}\left({\tilde{L}}_{O1}+{\tilde{L}}_{O2}\right)-\left({\tilde{C}}_{11}+{\tilde{C}}_{12}\right){\tilde{L}}_{O3}}{\left({\tilde{C}}_{11}-{\tilde{C}}_{12}\right)\left({\tilde{C}}_{11}+2{\tilde{C}}_{12}\right)},\\ {\tilde{C}}_{44}^{* } & = & {\tilde{C}}_{44}+2{\tilde{L}}_{23}{m}^{2}.\end{array}\end{eqnarray*}$$