Unified theory of magnetoelastic effects in B20 chiral magnets

A magnetic skyrmion is a spin whirl with topological protection and high mobility to electric current. Intrinsic magnetoelastic coupling in chiral magnets permits the manipulation of magnetic skyrmions and their lattice using mechanical loads, which is essential for developing future spintronics devices. It is found in experiments that the stability and deformation of skyrmions are sensitive to stresses, while the appearance of magnetic skyrmions in turn has a significant effect on the mechanical properties of the underlying material. However, a theory which explains these related phenomena within a unified framework is not seen. Here, we construct a thermodynamic model for B20 helimagnets incorporating a magnetoelastic functional with necessary higher-order interactions derived by group theory. Within the model, we establish the methodology to calculate the phase diagram and equilibrium properties of helimagnets under a coupled temperature-magnetoelastic field. Applying the model to bulk MnSi, we calculate the temperature-magnetic field phase diagram under stress-free condition and its variation when uniaxial compression is applied. We also calculate the variation of all the elastic constants with the magnetic field. The results obtained agree quantitatively with corresponding experiments. Our model provides a reliable basis for further theoretical studies concerning any magnetoelastic related phenomena in chiral magnets.


Introduction
Recent years have witnessed a revival of interest in chiral magnets such as MnSi, Fe 0.5 Co 0.5 Si, and FeGe, due to the experimental observation of a new chiral modulated magnetic state, commonly referred to as the skyrmion lattice phase [1][2][3]. The skyrmion lattice phase can be understood as the crystallization of isolated skyrmions, the latter of which are stabilized by the antisymmetric Dzyaloshinskii-Moriya (DM) interaction [4][5][6]. A skyrmion is attractive for its emergent electromagnetic properties such as spin motive force [7][8][9] and topological Hall effect [10,11]. The stability and appearance of skyrmions are sensitive to material size, shape [12][13][14][15][16] and various kinds of external fields [17][18][19]. Since the critical current density required to drive the motion of skyrmions is much lower than that for a magnetic domain wall [20,21], magnetic skyrmions are a promising candidate for the realization of the next generation spintronic devices.
Magnetoelastic coupling in chiral magnets permits interaction between skyrmions and the elastic fields of the underlying material, which leads to the occurrence of profound magnetoelastic phenomena. First, the application of mechanical loads can affect or even stabilize the skyrmion lattice. Through a theoretical model developed upon the Landau-Ginzburg functional, Butenko et al [22] find that uniaxial distortion stabilizes the skyrmion lattice in a broad range of thermodynamical parameters in cubic noncentrosymmetric ferromagnets, and further argue that this mechanism is responsible for the formation of skyrmion states observed in thin layers of Fe 0.5 Co 0.5 Si [2]. The stabilization of the skyrmion lattice is also observed in epitaxially grown FeGe thin films [23] and MnSi thin films [24], where the magnetoelastic effects should be significant due to the misfit strains exerted by the substrate. Later, it is observed in bulk MnSi that uniaxial tension affects the stability area of the skyrmion lattice phase in the temperature-magnetic field phase diagram [25]. Second, the existence of magnetic skyrmions in turn affects the elastic behavior of the material. Through ultrasonic studies [26,27], it is found that the appearance of the skyrmion lattice in chiral magnets is accompanied by a jump of elastic stiffness. Finally, deformation of the skyrmions and their lattice is strongly coupled with the deformation of the underlying materials. It is proved theoretically that the presence of skyrmions leads to nontrivial localized elastic fields [28,29] when the materials are free from external forces. And when uniaxial tension is applied, the skyrmion lattice in FeGe thin film [30] undergoes dramatic distortion two orders of magnitude larger than that of the underlying material.
To clarify the physical mechanism behind this, it is significant to establish a theory that can treat these different aspects of magnetoelastic phenomena within a unified framework. Previous theoretical studies on magnetoelastic coupling in MnSi or other chiral magnets mainly fall into two categories: one part is developed upon the magnetostriction theory constructed for ferromagnets with cubic symmetry (hereafter referred to as K theory) [31], while the other part is developed upon a Landau type mean-field model constructed for the specific spin-density-wave phase of MnSi (hereafter referred to as P theory) [32]. After thorough investigation, we find both theories to be oversimplified to explain skyrmion-related magnetoelastic phenomena; the variation of the elastic constants with the magnetic field examined for MnSi [26,27] in the skyrmion lattice phase cannot be understood within both theories. Recently, we became aware of a paper by Zhang and Nagaosa [33] addressing the ultrasonic elastic responses in a monopole lattice using an extended spin-wave theory concerning magnonphonon interaction. While such a microscopic model provides a deeper understanding of the origin of magnetoelastic coupling in chiral magnets, it is found that their model is more applicable to MnGe than MnSi. A possible reason is that the magnetoelastic Hamiltonian used in the model is constructed in a most simple form instead of a comprehensive description derived upon symmetry consideration. Besides, the effect of transverse acoustic waves is not taken into consideration in the model, with the result that it cannot be used to analyze the variation of the shear elastic constants C 44 and C 66 with the magnetic field.
In this paper, we formulate a thermodynamic model for B20 helimagnets incorporating a comprehensive magnetoelastic functional. The magnetoelastic functional is derived for the first time based on symmetry consideration of B20 helimagnets, where all necessary higher-order interactions are incorporated, so that variation of all the elastic constants of MnSi with the magnetic field observed in experiments can be quantitatively explained. The model is fundamental for studying any skyrmion-related magnetoelastic phenomena. Here, we explain two basic utilities of the model: (a) calculation of temperature-magnetic field phase diagram of any B20 compounds at any applied mechanical loads; (b) calculation of equilibrium properties for helimagnets at any given temperature, magnetic field, and mechanical loads. We apply the theory to bulk MnSi to evaluate the effect of magnetoelastic coupling on the temperature-magnetic field phase diagram, the equilibrium magnetization, and the elastic constants. By doing so, a complete set of thermodynamic parameters is determined for MnSi.
Neither the K theory nor the P theory is sufficient to explain the complex variation of elastic coefficients with the external magnetic field discovered in the experiments of MnSi [26,27], for which higher-order interactions w me1 andw me2 are introduced in equation (1).w me1 is needed when explaining the discrepancy between the elastic constants C 11 and C 33 in the skyrmion phase observed in an ultrasonic experiment of MnSi. On the other hand, w me2 is needed when explaining the variation of C 44 and C 66 with the magnetic field [27]. One should note that the two parts of the functional have already been simplified, where the details are described in appendix A.

Extended micromagnetic model incorporating magnetoelastic interaction
The Helmholtz free energy density for cubic helimagnets suffering coupled temperature-magnetoelastic field can be derived by incorporating the magnetoelastic functional developed in section 2 in the Ginzburg-Landau functional for chiral magnets with cubic symmetry [1,4,5,35,36]. We present the Helmholtz free energy density of the system in a rescaled form as denote, respectively, the rescaled anisotropic energy density and the rescaled elastic energy density, andw me is obtained by rescaling the magnetoelastic free energy density developed in section 2. In equation (7),w is given as a functional of the rescaled magnetization vector m and the elastic strains e ij at given rescaled temperature t and rescaled magnetic field b. Such a rescaled form reduces the number of effective thermodynamic parameters, and provides a material-independent theoretical framework to discuss the effect of magnetoelastic coupling. The rescaling process and the definition of all quantities with a wavy overline are given in appendix B. Equation (7) is fundamental to study various kinds of phenomena that occur in chiral magnets suffering coupled temperaturemagnetoelastic field. Here, we discuss two kinds of basic utilities as follows.
3.1. Temperature-magnetic field phase diagram calculation for helimagnets suffering mechanical loads To proceed, the first step is to solve the elastic strains at a given temperature, magnetic field, and mechanical loads, which differ for different kinds of mechanical boundary conditions. For the displacement boundary condition where the displacements are fixed at the boundaries as u , i0 the elastic strains are fully determined by the boundary condition as e e = ( ) u .  where s ij denotes the rescaled stress components, and * e ij and * g , ij the eigenstrains, are related to the rescaled magnetization by    In equations (11) and (12), the parameters with a superscript ' * ' are defined as * = Here, e s (˜) m, ij ij can be solved by taking the volume average of equation (10), while e ŝ (˜) m, ij ij can be derived by solving an eigenstrain problem [29]. For a mixed boundary condition, we generally have after deduction e e s = (˜) u m , , , where u i0 is the displacement prescribed at part of the boundary and s ij are the stresses solved using the stress boundary condition prescribed at the other part of the boundary. In all three cases, the elastic strains can generally be written as functions of the rescaled magnetization m: e e = ( ) m .

ij ij
In the second step, we substitute e e = ( ) m ij ij derived above into equation (7), which expresses the rescaled free energy density as a mere functional of m. Then, we consider a specific magnetic phase, which describes m by a certain mathematical expression. For chiral magnets, the known magnetic phases include: (i) the general conical phase  where q denotes the angle between the direction of q and [ ] 0 0 1 , sin sin cos . Here, it is pre-assumed that q always lies in the (110) plane, which is determined by the easy axis of the material [ ] 1 1 1 T and the direction of magnetic field [ ] 0 0 1 .
T It reduces to (ii) the conical phase when q = 0, which gives Equation (13) reduces to (iii) the general helical phase m ghelical when Equation (14) reduces to (iv) the helical phase m helical when = m 0, 3 and it reduces to (v) the ferromagnetic phase m ferro when = m 0. q (vi) The skyrmion lattice phase is described within the nth order Fourier representation [37] as: and n i denotes the number of reciprocal vectors whose modulus equals to s q.
i Here, s i is a positive sequence of numbers that increase with i, and m q ij can be expanded as (15) is constructed upon the hexagonal symmetry of the skyrmion crystal.
When distortion of the skyrmion crystal is considered, the emergent elastic strains have to be introduced in the Fourier representation [38]. By setting = n 1 and = = c c 0, 12 13 equation (15) reduces to the triple-Q representation which can be written without loss of generality as After specifying a magnetic phase, one solves the magnetization that minimizes the averaged free energy For example, if we consider the skyrmion lattice phase within the triple- , which determines the independent variables m m q , , In the third step, we repeat the free energy minimization process mentioned above for all possible magnetic phases. The equilibrium magnetic state is the one that yields the smallest averaged free energy density.
In the last step, we change the temperature and magnetic field, and then repeat the process mentioned above to determine the equilibrium magnetic state at the new condition. The temperature-magnetic field phase diagram is obtained after all points in the phase diagram are considered. A phase diagram of any two parameters can be derived in the same way if we change the temperature and magnetic field to two new parameters of interest.

Equilibrium properties for helimagnets concerning magnetoelastic coupling
To proceed, the first step is to determine the equilibrium magnetic state at a given temperature, magnetic field, and mechanical loads, using the method introduced above in part A. By doing so, we obtain the value of all independent variables at the given condition.
The second step is to decide which kind of equilibrium properties are to be discussed and clarify their thermodynamic definition, e.g. the rescaled elastic stiffness = = . Then we determine the work conjugates of all independent variables of the equilibrium magnetic state. The equilibrium property of interest is to be calculated at a given temperature and corresponding work conjugates of the independent variables. The third step is to solve the analytical expression of the equilibrium properties from the averaged free energy density by using the method of Jacobian transformation in thermodynamics [39]. The analytical expressions are generally lengthy and symbolic computation programs are needed for the deduction. When the analytical expressions are derived, we substitute the values of all independent variables obtained in the first step to calculate the values of the equilibrium properties.
Finally, we change the magnetic field and repeat the process above. The variation of the equilibrium properties of interest with the magnetic field is then obtained. Here, the magnetic field can be replaced by another parameter, and the variation of the equilibrium properties with the parameter can be obtained in the same way.
As an example, we derive for the first time the elastic stiffness in the skyrmion lattice phase at a given temperature, magnetic field and mechanical loads. At rescaled temperature t, rescaled magnetic field , , ,

Results for bulk MnSi
MnSi is a prototype material for us to understand the interaction between magnetic skyrmions and mechanical loads. In this section, we apply the theory established above to bulk MnSi. Within a unified theoretical framework, we are able to quantitatively explain four different aspects of the experimental results, including the temperature-magnetic field phase diagram for materials free from any mechanical loads, the variation of magnetostriction with magnetic field, the variation of elastic stiffness with magnetic field, and the temperaturemagnetic field phase diagram for materials suffering uniaxial pressure. By doing so, a comprehensive set of thermodynamic parameters for MnSi is obtained and listed in table 1, describing its magnetic, elastic, and magnetoelastic properties. The parameters are divided into two groups. The first group contains the parameters that can be directly obtained or have already been fitted from the experiments, while the second group contains the parameters that are fitted in this work. We briefly introduce the method and the experimental data used to fit the second group of parameters. From the magnetostriction experiment of bulk MnSi [28], L L L , , Oi can be fitted, where experimental results when the magnetic field is applied along three different directions (001), (110), and (111) are used. K L , 22 and L 23 are fitted from ultrasound measurements of the variation of the elastic coefficients with external magnetic field [27].
Using this set of thermodynamic parameters, extensive investigation is done concerning the magnetoelastic effects in MnSi, including the calculation of magnetostriction in the skyrmion phase and comparison to corresponding experiments, the calculation of the periodic elastic field in the skyrmion phase [29], the bumpy surface configuration of the skyrmion lattice [43], and the emergent elastic properties of the skyrmion lattice [38]. Here, we present three parts of the calculation results for bulk MnSi.   4.1. Temperature-magnetic field phase diagram concerning magnetoelastic coupling when the material is free from any mechanical loads Due to the magnetoelastic coupling, the elastic strains are related to the magnetization. Even when the system is free from any mechanical loads, we have nontrivial elastic strains from equation (10). This nontrivial e ij has an effect on the phase diagram as well as the solution of the equilibrium magnetization throughw el andw me in equation (7). In the phase diagram calculation of chiral magnets based on micromagnetic models [1,35,36], this back-action on the magnetization due to magnetoelastic coupling is usually neglected due to its smallness compared with other dominant terms in the free energy functional. To examine the applicability of such an assumption, we provide here a general analysis of the effect of magnetoelastic coupling on the magnetic phase diagram calculation for bulk chiral magnets when the system is free from any mechanical loads. We know that the coefficients of magnetoelastic coupling in different orders generally satisfy  ˜˜˜˜˜K L L L L q L q L q , , , , , 3 which gives * e e e µ K m , , 11 22 33 2 and that the influence of elastic strains on the equilibrium magnetization is predominantly attributed to the term ẽ Km ii 2 in equation (2). Hence, consideration of e e e , , 11 22 33 in the minimization of w renormalizes the coefficient of m 4 in the magnitude by -+˜. 11 12 In summary, we can compare the value of which suggests that magnetoelastic coupling should have a negligible effect on the shape of the magnetic phase diagram. For bulk MnSi, free from any mechanical loads, we plot the temperature-magnetic field phase diagram in figure 1(a) using the thermodynamic parameters fitted in table 1. The phase diagram is indistinguishable from the one plotted in our previous work [36] where magnetoelastic coupling is neglected. To further understand why this is the case, we plot in figure 1(b) the relative difference of free energy in the skyrmion phase - and the relative difference of free energy in the conical phase - Here, ( ) w m and ( ) w m 0 denote, respectively, the averaged free energy density calculated by integrating equation (7) and the averaged free energy density calculated by neglecting +w w el me in equation (7). We see that at 28 K, by including +w w el me in the free energy density, the free energy of the system in both the skyrmion phase and the conical phase decreases by about 2% to about 6% at different magnetic field. Yet the free energy in the two phases decreases in approximately the same way, which explains why the phase diagram is barely affected. On the other hand, magnetoelastic coupling pins the direction of the wave vectors (q , 11 q , 12 q 13 and so on) in the skyrmion lattice phase, as shown in figure 1(c). In the triple-Q representation of the skyrmion lattice, it is usually assumed that in equation (7)), the free energy density functional is invariant under an arbitrary rotation of q , 11 q 12 and q 13 in the x-y plane. After incorporating magnetoelastic coupling in the free energy density, we find that such a symmetry is broken and the direction of the triple-Q wave vectors is pinned. In figure 1(c), we plot in figure 2(a) the variation of  3  3  3  2  3  2 with j, where j describes the angle between q 11 and the x-axis (depicted in the inset of figure 2(a)). We find that at any temperature and magnetic field, the minimized free energy is found when j =  p . 2 As a result, the magnetoelastic coupling pins the triple-Q wave vectors to = [ ] q 0 1 0 ,   T m c c c  , , , , , , , , 3 is derived from equation (19), where ab C denotes a compressed matrix notation of the fourthorder tensor C ijkl [44]. The curves obtained from our theoretical calculation resemble corresponding experimental results in great detail [26,27]. To be more specific, in figure 2(a) DC 11 and DC 33 change to the opposite direction from the same point in the distorted conical phase as the magnetic field increases [26,27]; in the conical phase D > D C C 33 11 with a gap with a magnitude that approaches 10 GPa 8 and the gap mildly decreases as the magnetic field increases [26]; when a phase transition from the conical phase to the skyrmion phase occurs, we observe an obvious lift of DC 11 , which makes D > D C C 11 33 [26,27]; when a phase transition from the conical phase to the ferromagnetic phase occurs, DC 11 and DC 33 both increase, while DC 11 increases more sharply [27]. In figure 2(b), DC 44 and DC 66 change to the opposite direction in both the distorted conical phase and conical phase [27]; when a phase transition from the conical phase to the skyrmion phase occurs, DC 44 drops slightly, while DC 66 increases sharply [27]. In figure 2(c), variation of DC 12 and DC 13 with magnetic field is predicted, where corresponding experiments have never been performed before.
When plotting figure 2, the magnetization in the skyrmion phase is described by the third-order Fourier representation. We find that if we use the triple-Q representation instead, the values of D ab C obtained changes slightly within a range of 0.2%. This result shows that the order of Fourier representation of the skyrmion lattice phase has a negligible effect on the calculation of elastic constants. By substituting the solution of elastic strains into equation (7), we can initiate the where s ii is a stress tensor invariant, so that this result is valid for uniaxial compression applied in any direction. Hence, ẽ Km ii 2 renormalizes the second-order Landau expansion term in equation (7) as 11 12 for which the rescaled Curie temperature reduces from 1 to approximately s -+1 . K C C 2 11 12 For MnSi, uniaxial compression always decreases the Curie temperature, since <

Conclusion
In this paper, a thermodynamic model analyzing the coupled magnetoelastic fields in B20 helimagnets is developed based on group theoretical analysis. The model provides a unified theoretical framework to explain various aspects of skyrmion-related magnetoelastic experimental results, including but not limited to phase diagram calculation and equilibrium property calculation under coupled temperature-magnetoelastic field. By applying the model to bulk MnSi, we quantitatively reproduce the temperature-magnetic field phase diagram when the material is free from any mechanical loads, the variation of all the elastic constants with magnetic field, and the variation of temperature-magnetic field phase diagram when the material is suffering uniaxial compression in two different directions. We also obtain the general condition at which the effect of magnetoelastic coupling on the equilibrium properties can be neglected, and find that magnetoelastic coupling pins the triple-Q wave vectors of the skyrmion lattice phase in the x-y plane. Through calculation, we fit a whole set of thermodynamic parameters for MnSi, which lays a reliable foundation for further analytical or numerical analysis of magnetoelastic coupling phenomena. functional  Strictly speaking, point group T allows 13 independent thermodynamic parameters to describew .
T me1 Yet, in practice we do not have enough experimental data to fit all these parameters. Hence, we change the symmetry condition from point group T to point group O, which yieldsw me1 defined in equations (3) and (5). As forw , me2 we neglect all terms that are relevant to e , 11 e , 22 e , 33 becausew me2 describes higher-order effects of those described byw . me1 We find in calculation that these neglected terms relevant to e , 11 e , 22 e 33 are related to the variation of elastic constants such as C , 11 C . 33 But their contribution is negligible compared with corresponding terms inw .  (7), where