Thermodynamics of magnetic emergent crystals under coupled magnetoelastic fields

Magnetic fields and mechanical forces can change the deformation and stability of magnetic emergent crystals (MECs) such as Bloch skyrmion crystal (SkX), Néel SkX and Anti-SkX. Due to the tensor nature of strains, mechanical loads provide more fruitful ways to manipulate the MECs, while their effect on MECs other than the Bloch SkX is hitherto unclear. We construct a thermodynamic model for noncentrosymmetric ferromagnets in all possible point groups when subjected to coupled magnetoelastic fields. Compared with classic theories, we include terms coupling the elastic strains, the magnetization, and its derivatives in the free energy, which lead to strain-induced Dzyaloshinskii–Moriya interaction anisotropy. For epitaxial thin films in three types of point groups (T, C 3v , D 2d ) hosting Bloch SkX, Néel SkX and Anti-SkX, we find the newly added terms always deform the MECs and eventually lead to their instability as the misfit strains increase. Specifically, for Bloch SkX in group T materials and Néel SkX in group C 3v materials, a novel magnetic phase called paired-skyrmion crystal (pSkX) appears. Our theory lays the path to study deformation and phase transitions of different MECs, and to explore novel states of MECs in chiral magnets when subjected to magnetoelastic fields.


Introduction
Magnetic emergent crystals (MECs) are spatially periodic field patterns of magnetization that emerge from the atomic crystals. Examples of MECs include the Bloch skyrmion crystal (SkX) [1][2][3][4], Néel SkX [5][6][7], Anti-SkX [8,9], Meron crystals [10], etc. These magnetic states are attracting attention because of their ultrahigh mobility to electric current [11,12] which can be useful for racetrack memories [13][14][15], and because of the nano-scaled magnetic emergent particles, which permit bottom-up magnonics applications [16]. In both scenarios, the deformation of MECs is significantly affecting their desired properties [17][18][19][20][21]. When subjected to mechanical forces, SkX deforms significantly [22][23][24], affecting its generation [25], chirality [26], structural phase transitions [21], and elementary excitations [27]. The variety of related phenomena suggest that mechanical force is an effective approach to tune the properties of SkX and other MECs. Combining mechanical forces and magnetic fields, the tunability of the MECs is expected to be further enhanced. Practically, the presence of strains in materials is inevitable, either through possible defects and dislocations, or through lattice mismatch between the film and substrate in epitaxial thin films, or through surface tension in free nanomaterials. Therefore, a general theory which quantitatively determines the deformation of MECs and the variation of their physical properties when the material is subjected to magnetoelastic fields is in great need.
Theoretically, it is found that the magnetoelastic coupling of the underlying material provides a good explanation to the deformation of MECs in B20 compound [21,22], the phase diagrams [19], and variation of elastic constants of the materials [28,29] when subjected to mechanical forces [30]. In the model constructed, it is shown that classic theories of magnetoelastic coupling in ferromagnets [31,32] are not sufficient to explain the above mentioned phenomena, and terms of the form ε ij M k M l,n must be added to the Helmholtz free energy density. Interestingly, when studying the deformation of SkX in MnSi [21,30], it is found that the scale of the energy term ε ij M k M l,n is at least an order of magnitude smaller than that of the regular magnetoelastic term ε ij M 2 k , while the two terms are both indispensable. This recalls us of the fact that SkX in bulk helimagnets are stabilized by competition between the exchange interaction and the Dzyaloshinskii-Moriya interaction (DMI), where DMI is usually of a much smaller scale than that of the exchange interaction. In fact, the term ε ij M k M l,n can be interpreted as strain-induced DMI anisotropy, which generally exists for all chiral magnets, while their precise form for materials outside group T is not known. Besides, the free energy functional derived for B20 compounds is the Helmholtz free energy, which is applicable when the system is subjected to strain restrictions. When mechanical forces are applied at the boundary, the appropriate thermodynamic potential, being the elastic Gibbs energy [33], is not obtained.
In this paper, we systematically develop the thermodynamics of chiral magnets under coupled magnetoelastic fields. Specifically, we accomplish the following points: (A) for chiral magnets in different point groups, we derive the concrete form of magnetoelastic free energy density functional incorporating the strain-induced DMI anisotropy through symmetry analysis; (B) through Legendre transformation, we derive the general expression of the elastic Gibbs energy for chiral magnets; (C) for chiral magnets in three specific point groups (T, C 3v , D 2d ), we derive the thermodynamic potential for their epitaxial thin films when subjected to equiaxed in-plane misfit strain or pure shearing misfit strains. We systematically analyze the effect of the strain-induced DMI anisotropy on the deformation and stability of the three types of MECs (Bloch SkX, Néel SkX, Anti-SkX) that appear in the materials. Generally speaking, presence of the strain-induced DMI anisotropy term in the magnetoelastic free energy profoundly affects the deformation of the MECs such that they are stretched, sheared, or rotated in certain direction when subjected to a combination of magnetic fields and misfit strains. In particular, at appropriate conditions, a novel magnetic state called paired-skyrmion crystal (pSkX) may appear in thin films in group T and group C 3v . pSkX is also topologically nontrivial, but possesses a different topological density than the SkX phase. This means that as the value of misfit strains smoothly changes, a phase transition between two topologically nontrivial states occurs, accompanied by a jump of topological density. To sum up, we establish the theoretical foundation for studying magnetoelastic related phenomena in different skyrmionic magnetic materials.

Helmholtz free energy density with magnetoelastic interactions for chiral magnets
Generally, the Helmholtz free energy density with magnetoelastic interactions for chiral magnets takes the following form: where w mag (m) denotes the part of free energy density which relies merely on the magnetization vector m, w me (m, ε ij ) denotes the contribution of magnetoelastic coupling, and w el (ε ij ) denotes the elastic energy density which relies merely on the strain components ε ij . We present here two available models of w mag (m) that are widely used. The first model is temperature dependent, which yields a magnetization vector with changeable modulus: where w DM (m) denotes the generalized DMI whose exact form depends on the symmetry of the material (expanded in supplementary material (https://stacks.iop.org/NJP/23/023016/mmedia) part A), t denotes the temperature, and b denotes the magnetic field. This model derives from the Landau theory [32], which is usually used near the Curie temperature to study temperature dependent phenomena. The second model is used at low temperature, which yields a magnetization vector with fixed modulus: Equations (2) and (3) are presented in a rescaled manner for convenience, with two types of different rescaling procedure (for equation (2), see the supplementary information of reference [1], for equation (3), see reference [34]). For chiral magnets, we propose a general form of w me (m, ε ij ) as follow: where the first term with coefficient a ijkl (hereafter called the a term) is the interaction considered in classic theories of magnetoelastic coupling [31,32]. When studying the effect of strains in chiral magnets, the second term with coefficient d ijkln (hereafter called the d term) is non-ignorable, and sometimes dominant. If one fixes the strains of the material, the d term induces an anisotropy to the DMI, which can be significant due to the intrinsic smallness of the DMI. The form of magnetoelastic free energy density given in equation (4) can be simplified after considering the symmetry requirement of the materials, which reduces to where the first term with coefficient K Gi is reduced from the a term, and the second term with coefficient L Gi is reduced from the d term, and the subscript G depends on the point group of the material considered (e.g. for group O, the letter G is replaced by O). For a specific point group, n GK and n GL determines the independent thermodynamic parameters we have in the model. The detailed expression of f I Gi (ε ij , m) and f II Gi (ε ij , m, ∇m) are given in supplementary material part B. When the material studied is suffering homogeneous deformation, the number n GL can be further reduced to n GL , because some pairs of the expressions of f II Gi (ε ij , m, ∇m) are found to be equivalent to each other through integration by part. The values of n GK , n GL , and n GL are listed in table 1 for different groups.
To make use of the model proposed in equation (1), one has to obtain the exact values of the thermodynamic parameters, especially the magnetoelastic parameters a ijkl and d ijklm proposed in equation (4). These parameters can generally be determined through magnetostriction measurements. Here, we use group O as an example to show how to determine these thermodynamic parameters by experiments or first-principle calculation of the magnetostriction. For materials in group O, we have where When the material studied is free from external forces, we have for the stresses σ ij = ∂(w me +w el ) 44 (ε 2 23 + ε 2 13 + ε 2 12 ) + C 12 (ε 11 ε 22 + ε 11 ε 33 + ε 22 ε 33 ). For materials in group O it gives: where C ij are the elastic stiffness of the material, and To determine the parameters, one should apply bias magnetic field of different strength in different directions, and then measure the associated strains ε ij . Firstly, apply a strong magnetic field in axis-r 3 at low temperature, the magnetization in such a case should be m = [0 0 m s ] T , where m s denotes the saturation magnetization. For the measured strains ε ij , we have after manipulation Equations (9) and (10) Similarly, apply a strong magnetic field and a moderate magnetic field in direction [110], the magnetization is given respectively by We have after manipulation Equations (9)-(15) provide the relations between the thermodynamic parameters and the measured magnetostriction (in equation (11), L O1 − L O2 is treated as a single parameter because when suffering homogeneous elastic fields, the terms related to L O1 and −L O2 can be switched to each other through integration by part; equation (15) is understood in the same way). Further study shows that the term with coefficient L O4 in equation (6) is not affecting the magnetostriction measurement. It only has an effect on the inhomogeneous distribution of periodic strains in space, for which the value of parameter L O4 cannot be examined through ordinary magnetostriction measurement. For materials in points groups with lower symmetry, we need to perform magnetostriction measurements in more directions to get a full set of parameters and the solution of the parameters is generally more complicated.

Elastic Gibbs energy density with magnetoelastic interactions for chiral magnets
In the Helmholtz free energy density given in equation (1), the elastic strains and the magnetization are chosen as the independent variables of the system, which is applicable to materials with a fixed displacement, or approximately, with displacement boundary conditions. When the stress boundary conditions are given instead, the elastic Gibbs free energy density [33] is the appropriate thermodynamic potential to be used, which can be derived from the Helmholtz free energy density through the Legendre transformation. We derive it as follow for chiral magnets suffering a coupled magnetoelastic field. From equation (4), all strain-related terms in w me are linear to the strain components. Therefore equation (4) can be reformulated as where The elastic Gibbs energy density of the system can be derived through the following Legendre transformation where σ = [σ 11 σ 22 σ 33 σ 23 σ 13 σ 12 ] T . After manipulation, we have from equation (19) g where w el (σ) = 1 2 σ T Sσ, and S = C −1 is the flexibility matrix. Equation (20) is a general formula of the elastic Gibbs energy density. To use it, one has to first consider the specific symmetry requirement of the material studied, and derive the particular form of F from equation (5), and then substitute it into equation (20). Notice that the last term on the right-hand-side of equation (20) is generally very complicated, while its effect is to induce a modification of the magneto-crystalline anisotropy and other higher order magnetic anisotropy terms. Yet according to previous analysis [30], for materials with a magnetostriction smaller than 10 −4 (which is generally true), the effect of this term is negligible. For different point groups, the detailed expressions of F ij (m, ∇m) are given in supplementary material part C.

Thermodynamic potentials for epitaxial thin films of chiral magnets
Sometimes the material studied possess a mixed boundary condition. The most frequently encountered case of this type is an epitaxial thin film, which is mechanically clamped in-plane, and subjected to forces on the out-of-plane surfaces. In this case, the appropriate thermodynamic potential is: where the independent variable for the elastic fields are ε 11 , ε 22 , ε 12 , σ 33 , σ 23 and σ 13 . The expression of g TF depends on the symmetry of the material studied. For orthorhombic materials (group D 2 or C 2v ), we have where g el denotes the part of free energy density that depends merely on the elastic strains, and in the case of an epitaxial strain is equal to a constant, for which its detailed expression is not expanded. The expressions of F ij differ for different groups and can be found in supplementary material part C.
For trigonal materials (group D 3 , C 3 or C 3v ), we have For hexagonal materials (group D 6 , C 6 or C 6v ), we have For cubic materials (group O or T), we have
For materials in group T subjected to condition (a), the thermodynamic potential to be used is reduced to , and K Ti and L Ti are the magnetoelastic thermodynamic parameters for materials in group T, which are defined in supplementary material part B. When K * T1 = 0 and L * T2 = −L * T1 , equation (26) reduces to the thermodynamic potential for materials in group O. Equation (26) shows that for epitaxial helimagnets in group T subjected to equiaxed in-plane misfit strain, the emergent deformation of the MECs are completely determined by six thermodynamic parameters: K * T , K * T1 , K * T2 , L * T1 , L * T2 , and L * T3 . Here we investigate the effect of each individual magnetoelastic terms in equation (26). K * T εm 2 reduces the Curie temperature of the material by K * T ε [here and in the discussion below we have used equation (2) instead of equation (3)]. K * T1 εm 2 1 and K * T2 εm 2 3 induce an uniaxial anisotropy, which is studied in previous works [41,42]. The term with a coefficient L * T1 is ineffective when studying 2D MECs in the r 1 -r 2 plane. The term with a coefficient L * T2 shares the form of the original DMI of the material, for which its effect is to change the lattice constant of the MEC studied without causing a pattern deformation. The term with a coefficient L * T3 induces an anisotropy of DMI, which causes deformation and even phase transition of MECs (figure 1). To be more specific, a negative (positive) L * T3 ε stretches the Bloch SkX in direction r 1 (r 2 ), and as the value of L * T3 ε decreases (increases), the Bloch SkX finally evolves to the In-plane single Q (IPSQ) phase.
One should notice that in a certain range of magnetic field, as the value of positive L * T3 ε increases, a 'cell division' phenomena occurs: the skyrmions in the SkX is suddenly stretched into a pair of 'linked' skyrmions, which form a novel magnetic state called pSkX. The 'particle' of pSkX, a paired-skyrmion, is different from a biskyrmion [43] or a dipole skyrmion [44]. To be more specific, a paired-skyrmion is composed of two identical skyrmions; a biskyrmion is composed of two skyrmions with the opposite in-plane magnetization distribution and the same out-of-plane magnetization distribution; a dipole skyrmion is composed of two skyrmions with the opposite in-plane and out-of-plane magnetization distribution. Interestingly, the topological density of pSkX is nontrivial, but is different from that of the SkX state. As a result, the topological density of the system undergoes a sudden change during a SkX-pSkX phase transition, which is realized by smoothly changing the applied misfit strains. The details for calculating the phase diagrams and magnetization distributions in figure 1 and other figures are introduced in the methods section.
For materials in group T subjected to condition (b), the thermodynamic potential to be used is reduced to g TF = w mag + g el + 2K T4 γm 1 m 2 + L T7 γm 1 m 1,3 + L T9 γm 2 m 2,3 + L T13 γm 3 m 3,3 where L * T4 = L T12 − L T10 , L * T5 = L T11 − L T8 . When L T9 = −L T7 , L T13 = 0, and L T11 − L T8 = L T10 − L T12 , equation (27) reduces to the thermodynamic potential for materials in group O. Here we investigate the effect of each individual terms in equation (27). The term with a coefficient K T4 can be transformed to uniaxial anisotropy in a particular direction in the r 1 -r 2 plane (i.e. after a rotation of coordinates in the r 1 -r 2 plane, the term with a coefficient K T4 can take a similar form to K * T1 εm 2 1 in equation (26), for which they possess similar effects). Generally speaking, for all the point groups studied, all terms with a coefficient K Gi can be transformed to uniaxial anisotropy in a specific direction, for which we will not analyze the effect of terms with a coefficient K Gi repeatedly in the discussion for materials in other point groups. Terms with coefficients L T7 , L T9 , and L T13 are ineffective when studying 2D MECs in the r 1 -r 2 plane. The two terms with coefficients L * T4 and L * T5 induce anisotropy of DMI, which cause deformation and even phase transition of MECs (figure 1). To be more specific, a negative (positive) L * T4 γ stretches the Bloch SkX in direction [110] ([110]), and as the value of L * T4 γ decreases (increases), the Bloch SkX finally evolves to the IPSQ phase. The term with a coefficient L * T5 has an exactly opposite effect to that of the term with a coefficient L * T4 (if we assume that L * T5 = −L * T4 , the L * T5 term has exactly the same effect as the L * T4 term). For materials in group C 3v subjected to condition (a), the thermodynamic potential to be used is reduced to g TF = w mag + g el + K * C 3v εm 2 + K * C 3v 1 εm 2 3 + L * C 3v 1 ε(m 1 m 1,1 − m 2 m 2,1 ) + L * C 3v 2 ε(m 1 m 1,3 + m 2 m 2,3 ) where , L * C 3v 4 = L C 3v 10 + L C 3v 12 − L C 3v 4 − L C 3v 9 + 2C 13 C 33 L C 3v 27 − 2C 13 C 33 L C 3v 26 , and K C 3v i and L C 3v i are the magnetoelastic thermodynamic parameters for materials in group C 3v , which are defined in supplementary material part B. If L * C 3v 1 = 0, equation (28) reduces to the thermodynamic potential for materials in group C 6v . Here we investigate the effect of each individual terms in equation (28). Terms with coefficients L * C 3v 1 , L * C 3v 2 , and L * C 3v 3 are ineffective when studying 2D MECs in the r 1 -r 2 plane. The term with a coefficient L * C 3v 4 shares the form of the original DMI of the material, for which its effect is to change the lattice constant of the MEC studied without causing a pattern deformation.
For materials in group C 3v subjected to condition (b), the thermodynamic potential to be used is reduced to where g n = L * C 3v 9 γm 1 m 1,2 + L * C 3v 10 γm 2 m 2,2 − and , L * C 3v 9 = 1 4 (−L C 3v 1 + L C 3v 3 + L C 3v 6 − 3L C 3v 7 ), L * C 6v 10 = − 1 4 (3L C 3v 1 + L C 3v 3 + L C 3v 6 + L C 3v 7 ). When K C 3v 2 = K C 3v 5 = K C 3v 6 = L * C 3v 1 = L * C 3v 5 = L * C 3v 6 = L * C 3v 8 = 0 and g n = 0, equation (29) reduces to the thermodynamic potential for materials in group C 6v . Terms in equation (30) are ineffective when studying 2D MECs in the r 1 -r 2 plane. The four terms with coefficients L * C 3v 5 , L * C 3v 6 , L * C 3v 7 and L * C 3v 8 induce anisotropy of DMI, which cause deformation and even phase transition of MECs (figure 2). To be more specific, as the modulus of L * C 3v 5 γ increases, the Néel SkX rotates about 30 • . The term with a coefficient L * C 3v 6 causes internal deformation [45] of the Néel SkX: as the modulus of L * C 3v 6 γ increases, the Néel skyrmions in SkX are stretched in axis r 1 , while the lattice constant of the Néel SkX is kept unchanged. The term with a coefficient L * C 3v 7 causes pure shearing of the Néel SkX: a positive (negative) L * C 3v 7 γ rotates and stretches the SkX in direction [110] ([110]). The term with a coefficient L * C 3v 8 causes uniaxial elongation of the SkX: a positive (negative) L * C 3v 8 γ stretches the SkX in axis r 1 (r 2 ). One should notice that in a certain range of magnetic field, as the value of negative L * C 3v 8 γ decreases, the SkX-pSkX phase transition with is observed in Bloch SkX for materials in group T occurs again. Similarly, the topological density of the system undergoes a sudden change during a SkX-pSkX phase transition of materials in group C 3v .
For materials in group D 2d subjected to condition (a), the thermodynamic potential to be used is reduced to C 33 L D 2d 20 , and K D 2d i and L D 2d i are the magnetoelastic thermodynamic parameters for materials in group D 2d . In equation (31), the term with a coefficient L * D 2d 1 is ineffective when studying 2D MECs in the r 1 -r 2 plane. The term with a coefficient L * D 2d 2 shares the form of the original DMI of the material, for which its effect is to change the lattice constant of the MEC studied without causing a pattern deformation.
For materials in group D 2d subjected to condition (b), the thermodynamic potential to be used is reduced to (e) plots the variation of the topological density of magnetization field with b corresponding to the loading condition of (d), and the result is calculated at λ C 3v 4 = −0.2. In (a i , b i , c i , d j ), the vectors illustrate the distribution of the in-plane magnetization components with length proportional to their magnitude, while the colored density plot illustrates the distribution of the out-of-plane magnetization component. The black where L * D 2d 3 = (L D 2d 9 − L D 2d 8 ). In equation (32), terms with coefficients L D 2d 7 and L D 2d 10 are ineffective when studying 2D MECs in the r 1 -r 2 plane. The term with a coefficient L * D 2d 3 causes pure shearing of the Anti-SkX  To sum up, we systematically analyze the effect of each of the strain-induced DMI terms on the deformation and phase transition of MECs appearing in chiral magnets in group T, C 3v , and D 2d when the thin film is subjected to misfit strains of equiaxed in-plane tension/compression or pure shear. In real materials, the condition of misfit strains is usually more complicated, and all the strain-induced DMI terms coexist with different strength depending on the material studied. The analysis performed in this work can be regarded as basic case studies of the large variety of possible deformation pattern and novel MECs that can appear due to presence of coupled magnetoelastic fields.

Methods
In this section, we introduce how to use the free energy density derived in the above section to perform free energy minimization for different MECs. In particular, we discuss the terms in the magnetoelastic free energy density which may lead to anisotropic deformation or phase transition of the MECs. For epitaxial thin films in group T subjected to condition (a), we should discuss the equilibrium state of the material by using equation (26). Based on the discussion given below equation (26), the only magnetoelastic term that leads to anisotropic deformation of the Bloch-SkX is the term with coefficient L * T3 . To study the effect of this L * T3 term alone, we assume that λ T1 = L * T3 ε, and K * T = K * T1 = K * T2 = L * T1 = L * T2 = 0, t = 0.5 in equation (26), for different values of λ T1 and the magnetic field b = [00b] T , we minimize the volume integral of equation (26) to determine the equilibrium magnetic state. And by changing the values of λ T1 and b, we repeat the minimization process which finally yields the λ T1 -b phase diagram as shown in figure 1(a). Here the variation of the parameter λ T1 has two ways of interpretation: it can either be understood as the change of misfit strain ε for a certain type of material, or be understood as a variation of the parameter L * T3 while the misfit strain is fixed. When practically performing the calculation, we choose w mag in equation (26) to be w mag1 (m) defined in equation (2). We consider the case where the thickness of the thin film is smaller than a period of the helical state of the material, in which case the conical state is suppressed, so that the magnetization of the material can be expression by the Fourier representation as [45,46]: where q ij denotes the reciprocal vectors of the undeformed EC, which is ordered in such a way that q 1j < q 2j < q 3j < · · · , and |q i1 | = |q i2 | = · · · = q in i . In equation (33), m 0 is a constant vector, and F e = ε e 11 ε e 12 + ω e ε e 12 − ω e ε e 22 (34) with ε e ij being the emergent elastic strains of the EC, and ω e being the emergent rotational angle. In practice, we truncate the expansion given in equation (33) at n = 3, and substitute it into equation (26) (one should notice that when n = 3, we have considered 18 Fourier terms in equation (33) while the well-known triple-Q approximation considers only 6 terms). At given values of λ T1 and b, we minimize the free energy with respect to m 0 , m q ij , ε e ij , and ω e , which determines the equilibrium magnetization. Repeating this free energy minimization process for different values of λ T1 and b, one obtains the variation of magnetization with them, as well as the corresponding phase diagrams.
The calculation process introduced above is generally useful. For epitaxial thin film of materials in group T subjected to condition (b), we consider the following thermodynamic parameters: λ T2 = L * T4 γ, λ T3 = L * T5 γ, and b, and minimize the free energy given in equation (27). For epitaxial thin film of materials in group C 3v subjected to condition (b), we consider the following thermodynamic parameters: λ C 3v 1 = L * C 3v 5 γ, λ C 3v 2 = L * C 3v 6 γ, λ C 3v 3 = L * C 3v 7 γ, λ C 3v 4 = L * C 3v 8 γ and b, and minimize the free energy given in equation (29). For epitaxial thin film of materials in group D 2d subjected to condition (b), we consider the following thermodynamic parameters: λ D 2d 1 = L * D 2d 3 γ and b, and minimize the free energy given in equation (31). Softwares that permit symbolic calculations are helpful when doing the derivation, and we have used mathematica in this work.

Data availability statement
All data that support the findings of this study are included within the article (and any supplementary files).