Wave-nature and metastability of emergent crystals in chiral magnets

Topological spin textures emerging in magnetic materials usually appear in crystalline states. A long-standing dilemma is whether we should understand these emergent crystals as gathering particles or coupling waves, the answer of which affects almost every aspect of our understanding on the subject. Here we prove that 2-D emergent crystals with long-range order in helimagnets, such as skyrmion crystals and dipole skyrmion crystals, have a wave nature. We systematically study their equilibrium properties, metastability, and phase transition path when unstable. We show that the robustness of a skyrmion crystal derives from its metastability, and that its phase transition dynamics at low (high) magnetic field is mediated by a soft mode which breaks (maintains) its hexagonal symmetry. Different from ordinary crystals which are formed by. and breaks into atoms, emergent crystals have a new formation (destruction) mechanism: they appear from (turn to) single-Q spin-density-wave states through nonlinear mode-mode interactions.

M agnetic skyrmions 1-3 are topological spin textures emerging in magnetic materials due to non-collinear interactions permitted by inversion symmetry breaking. They are attractive due to their high mobility to electric current 4,5 , tunability to material size 6,7 , and sensitivity to various kinds of field [8][9][10][11][12] . Besides ordinary skyrmions, other magnetic "emergent particles" have also been discovered in recent years: anti-skyrmions 13 are found in non-centrosymmetric ferromagnets (antiferromagnets) with D 2d and S 4 symmetry, biskyrmions 14 are stabilized by demagnetization effects in common magnetic materials, and dipole-skyrmions 15,16 are found in magnetic multilayers with dipole interactions. Instead of isolated individuals, these emergent particles generally appear in crystalline states with long-range order 17 in most bulk materials 1,10,18 and layered structures [14][15][16][19][20][21][22][23][24][25] where they are detected. A fundamental question under debate is: should we understand emergent crystals as gathering particles [26][27][28][29] or coupling waves 1,17,18 ? Since the local symmetry permitted by the field configuration of an emergent particle is usually incompatible with the global symmetry permitted by its emergent crystalline state, the answer to this question is vital for us to understand the basic properties of emergent crystals and their relation with the composing emergent particles.
In this work, we prove based on strict variational analysis that 2-D emergent crystals appearing in helimagnets with long-range order have a wave nature. Based on a wave description of four types of emergent crystals, we not only obtain their equilibrium properties, but also successfully reproduce their phase transition dynamics at both low and high magnetic fields. The phase transition behavior of skyrmion crystal (SkX) at low magnetic field reveals a novel formation (destruction) mechanism of emergent crystals, such that they emerge from (turn to) single-Q spindensity-wave states through non linear mode-mode interactions. This mechanism is not seen for ordinary crystals.

Results
Variational analysis of emergent crystals in helimagnets. We perform our analysis based on a rescaled free energy density functional 1,27,30,31 for cubic helimagnets (see the Methods section for details)w where m, t, and b denote respectively the rescaled magnetization vector, the rescaled temperature and the rescaled magnetic field.
An emergent crystal appears as a periodic solution of m(r i ) that minimizes the averaged free energy density ðmÞdV. At first glance, it is straightforward to expand m(r i ) in a Fourier series if we presuppose the existence of a periodic magnetization structure. However, when a function has discontinuous first order partial derivatives at certain points, the "Gibbs phenomenon" 32 is induced which makes it not suitable for Fourier analysis (see a detailed Fourier analysis of the circular cell approximation (CCA) 26 in Supplementary Note 1 and 2, where one should keep in mind that the CCA configuration keeps simultaneously the local axial symmetry of an isolated skyrmion and the global hexagonal symmetry of the lattice).
In this case, a proper question to ask is: to be a solution that minimizes wðmÞ, can m have discontinuous first order partial derivatives at certain boundary in the unit cell? Concerning 2-D emergent particles only, the shape of this boundary should be a circle so that inside the circle the axial symmetry of m is maintained due to the axial symmetry of Eq. (1), and outside the circle the axial symmetry of m is destroyed due to long-range order of the emergent crystal. We construct the following field configuration in Fig. 1: in a hexagonal unit cell, m is axially symmetric for ρ ≤ r 0 (in the orange region), and is not axial symmetric for ρ > r 0 (in the blue region). The integrability condition of wðmÞ requires that m is continuous everywhere, and in both regions 0 < ρ < r 0 and ρ > r 0 , the first order derivatives of m with respect to the spatial variables are continuous. It is convenient to describe m in spherical coordinates by (m, θ, ψ) and to describe the space in cylindrical coordinates (ρ, φ, z). At ρ = 0, local axial symmetry requires sinθ = 0, while ψ ρ ¼ ∂ψ ∂ρ does not need to be continuous. Since m is continuous on the interface ρ = r 0 , m φ , θ φ , and ψ φ are all continuous on the interface ρ = r 0 . Such a field configuration of m keeps simultaneously the local axial symmetry in ρ ≤ r 0 and the global hexagonal symmetry. Yet, for such a field configuration of m to be the solution of the variational problem, additional conditions are required. Applying the Haar's theorem 32 in an arbitrary sector region K depicted in Fig. 1 to the variational problem that minimizes wðmÞ, we find that the solution function m = m(ρ, φ), θ = θ(ρ, φ), and ψ = ψ(ρ, φ) must satisfy where λ is to be replaced by m, θ, or ψ. Let r 1 = r 0 − ε, r 2 = r 0 + ε, and substituting Eq. (1) into Eq. (2), for ε → 0 we have after manipulation Fig. 1 Illustration of a hexagonal unit cell Ω. In the orange circular region bounded by ρ = r 0 , the field configuration of m is axially symmetric, while in the rest of the cell (blue region), the axial symmetry is broken. We construct a sector region K bounded by ρ = r 1 , ρ = r 2 , φ = α, and φ = β where we have used the fact that α and β are arbitrary. Equations (3) and (4) show that m ρ and θ ρ are continuous on the interface ρ = r 0 . Things are more complicated for Eq. (5), where we have at the interface ρ = r 0 either sinθ = 0 or ψ ρ to be continuous. Nevertheless, we find after derivation that in both cases ((1) sinθ = 0 and ψ ρ incontinuous; (2) ψ ρ continuous), ∇m and ∇ × m are continuous at ρ = 0 and ρ = r 0 , which makes them continuous everywhere. For a continuous periodic vector function v, its Fourier series expansion v Fn is uniformly convergent to v as n → ∞. As a result, the solution function m min of the variational problem, as well as ∇m min and ∇ × m min can all be expressed as Fourier series. Since Fourier series expansion and differentiation commute 33 , they are uniformly approached by m mFn , ∇m mFn , and ∇ × m mFn , where m mFn is the Fourier series expansion of m min . The minimized rescaled free energy densitywðm min Þ, as a continuous function of m min , ∇m min , and ∇ × m min , is therefore uniformly approached bywðm mFn Þ. Thus the Fourier representation is the exact form of solution to describe any emergent crystals in chiral magnets.
Fourier representation of an emergent crystalline state. The proof given above leads naturally to the following two points: (a) crystallization of isolated skyrmions into SkX is a phase transition, characterized by breaking of the local axial symmetry of an isolated skyrmion. (b) In general, the order parameter of any emergent crystals (e.g., m of SkX) can be described by the nth order Fourier representation as where jq i1 j ¼ jq i2 j ¼ jq i3 j ¼ ¼ jq in i j ¼ s i q; jq 1j j <jq 2j <jq 3j j < <jq nj j, and n i denotes the number of reciprocal vectors whose modulus equal to s i q. One should notice that such a description is generally applicable to any n-dimensional emergent crystals (n = 1, 2, 3). Some specific information about the Fourier representation of 2-D hexagonal lattices is listed inTable 1. We describe four types of emergent crystals here using Eq. Specific attention should be paid on the modulus of magnetization given in Eq. (6). Within a truncated Fourier representation, we find that the modulus of magnetization can never be a constant in space if m 0 and m q ij are all non-zero vectors. For example, the term P n i j¼1 m 0 Á m q 1j e iq 1j Ár appearing in the expansion of |m Fn | 2 can never be canceled by increasing the Fourier expansion order. Further analysis tells us that |m Fn | 2 approaches a constant only when all m q ij approaches zero (ferromagnetic phase). This modulations of the magnetization modulus is found to be important to understand the existence of precursor states near the ordering temperature discovered in previous studies within the functional given in Eq. (1) 34,35 .
Equilibrium properties of emergent crystals. The Fourier representation provides an analytical framework to study the equilibrium properties of any emergent crystal. For a specific type of emergent crystal, one calculates the partial derivative of wðmÞ with respect to m 0 , m q ij , and q, which yields a set of high order algebraic equations. For the three types of non-trivial emergent crystals, we calculate the variation of wðmÞ with b at t = 0.5 and compare it with that for the conical phase in Fig. 3a. Within mean-field approximation, the conical phase always has the lowest averaged free energy density, while DiskX always has the highest averaged free energy density. For all emergent crystals, the averaged free energy density is always slightly decreased for an increase of the Fourier representation order, in which case 1st order Fourier representation is found to be an effective approximation. On the other hand, the variation of lattice constant of emergent crystals with b is found to be very sensitive to Fourier representation order (Fig. 3b). To be more specific, the m 4 term in Eq. (1) permits non-zero mode-mode interactions, such as that described by m 0 Á m q 3k m q 1j Á m q 1j , which is negative for appropriate choice of m q 3k . Meanwhile, non-zero m q 3k reduces the equilibrium value of q through the terms P 3 i¼1 ∂m ∂r i 2 þ2m Á ð∇ mÞ in Eq. (1). As the reciprocal of q, the lattice constant generally increases with any m q ij ; ði ! 2Þ. For n ≥ 3, the variation of lattice constant of SkX with magnetic field obtained in experiment 33 is well explained within the nth order Fourier representation of the SkX. The equilibrium properties of the SkX phase is determined through numerically solving a set of algebraic equations, where the equations are derived analytically based on the Fourier representation of SkX. Nevertheless, this Fourier-based method is fundamentally different from pure numerical methods for free energy minimization such as the finite difference method and the Monte Carlo simulation method, since it provides the analytical expression of the magnetization. Due to the lack of an analytical expression in pure numerical methods, all information of SkX has to be extracted vaguely from the numerical distribution of magnetization, which is sometimes unreliable. For instance, the core region of a skyrmion in the SkX phase was thought to possess axial symmetry according to finite difference calculation 35 . From the Fourier representation of SkX, we learn that this is incorrect once SkX maintains strict long-range order, because a Fourier series destroys axial symmetry regardless of the expansion order. Moreover, similar to pure numerical methods, the precision of the Fourier-based method can be improved by increasing the Fourier expansion order, and when achieving the same level of precision, we find that the Fourier-based method is much more efficient than the Monte Carlo program we developed (See Supplementary Note 3 for more details).
Metastability of emergent crystals. One exotic property of SkX is that it is sometimes found to be stable at conditions where it is not supposed to be 18 . We explain this robustness of SkX by its metastability: SkX corresponds to a local minimum of the free energy surrounded by high barriers, and the thermal fluctuation is not strong enough for the system to escape from this trap. In this sense, the critical condition when emergent crystals losses their metastability provides an upper bound concerning their robustness. In the Methods section, we establish the method of soft-mode analysis to determine this critical condition for emergent crystals. Applying the method to SkX at the condition t = 0.7, b > 0, we find that the metastability is dominated by two different modes at low and high magnetic field, which turn to soft modes at two critical magnetic fields b csl and b csh (Fig. 4a, b). The variation of b csl and b csh with t determines the region of metastability of SkX in the t-b phase diagram (Fig. 4c). Similar analysis to the IPSQ phase and ferromagnetic phase yields two different critical magnetic fields b ci and b cf . The discrepancy between b ci and b csl (b cf and b csh ) provides a region in the t-b phase diagram (shadowed area in Fig. 4c) where both the IPSQ phase and SkX (the ferromagnet phase and SkX) are metastable, so no phase transition between the two states will spontaneously occur, and a two-phase state may be observed. Such a two-phase state of the IPSQ phase and SkX (the ferromagnet phase and SkX) is observed in TEM experiments of Fe 0.5 Co 0.5 Si 19 and FeGe 21 thin films, while a two-phase state of the DiskX and the IPSQ phase is also observed in La 2 −2x Sr 1 +2x Mn 2 O 7 14 and (Mn 1−x Ni x ) 65 Ga 35 25 . In a two-phase state, SkX and the IPSQ phase occupy different parts of the space, forming domain-like structures. Near the domain walls, intermediate state such as those depicted in Fig. 2f will appear locally.
It is shown in Fig. 3b that the most stable phase always corresponds to the conical phase at the condition t = 0.5,0 < b < 0.25. We find that the conical phase is approached through a spin reorientation transition from the IPSQ phase. To account for this transition, one has to introduce the wave-vector-field angle into the expression of single-Q magnetization to describe the so-called general conical phase 31 . One should notice that although the conical phase corresponds to the state with global minimum of free energy at low magnetic field, the SkX phase always evolves to the IPSQ phase first according to the soft-mode analysis. This is because in the free energy landscape, the IPSQ phase is "closer" to the SkX phase than the conical phase in the sense that the first two states share the same form of Fourier representation.
Phase transition dynamics of skyrmion crystal. Another interesting phenomenon about SkX is the distinct phase transition behavior it shows at low and high magnetic field. TEM experiments of Fe 0.5 Co 0.5 Si 19 and FeGe 21 thin films both show that SkX becomes elongated in a certain direction at enough low magnetic field and finally evolves to the IPSQ phase, while at high enough magnetic field it evolves to the ferromagnetic phase without a distortion of the skyrmions. We reproduce the two phase transition paths of SkX explicitly in Fig. 2. The two different transition paths can be understood by investigating the characteristic vectors V C of the corresponding soft modes. According to its definition, a phase transition of SkX maintains the hexagonal symmetry if V C of the corresponding soft mode satisfies (V C ) 1 =  6 , which explains the different phase transition paths of SkX at low and high magnetic fields shown in Fig. 2. Specifically, the mechanism for the SkX-IPSQ transition lies in the non-linear mode-mode interaction. Within the 1 st order Fourier representation, the triple-Q SkX induces two more types of coupling compared with the IPSQ phase due to the m 4 term in the free energy density, one represented by m 2 q 11 m 2 q 12 , and the other represented by m 0 Á m q 13 m q 11 Á m q 12 . The first type of interaction always increases the free energy density of SkX compared with that of the IPSQ phase, while the second type of interaction can effectively reduce the free energy density of SkX by appropriately choosing the direction between m 0 and m q 11 ; m q 12 ; m q 13 . Since |m 0 | → 0 when b → 0, the second type of interaction has a negligible effect when the magnetic field is low enough, for which the IPSQ phase is more stable than SkX. As b increases, the second type of interaction becomes more significant, and at certain point SkX becomes more stable. Since interaction between waves occurs globally, the stability of SkX at low magnetic field is found to be irrelevant to the local topological protection of an isolated skyrmion. Moreover, it implies that emergent crystals composed of emergent particles that are unstable alone may appear under appropriate conditions, which is supported by the phase transition from a triangular SkX to a square SkX observed in a Co-Z-Mn alloy 18 .
Applying the method of soft-mode analysis to DiskX, we find that DiskX is always unstable in chiral magnets (Fig. 4) and spontaneously evolve to SkX with a phase transition path shown in Fig. 2a, c, d. The triangle chiral whirl with magnetization pointing downwards (green whirl in Fig. 2a) will slightly enlarge and finally become a skyrmion with hexagonal shape in SkX.
Wave-particle duality for emergent crystals. The analysis above shows that the wave description of an emergent crystal is not only mathematically correct, but also physically indispensable. The phase transition of SkX to the IPSQ state at low magnetic field, mediated by softening of macroscopic spin-density-waves, cannot be understood by merely studying the localized skyrmion-skyrmion interaction. Moreover, the wave description of an emergent crystal does not deny the particle nature of its localized composing units. The phase transition of SkX to the ferromagnetic phase at high magnetic field is mediated by a diverging expansion of the emergent crystal lattice which maintains the crystalline symmetry (Fig. 2g). This result coincides with previous analysis 34, 35 , indicating a crossover phenomenon of the skyrmion-skyrmion interaction from mutual attractive to mutual repulsive. In this process, isolated skyrmion should appear as an intermediate state when the skyrmion-skyrmion interaction is weak enough. These two phase transitions show precisely the wave-particle duality of the SkX phase, which is also a fundamental difference between emergent crystals and atomic crystals. Unlike atomic crystals which are always condensates of atoms, emergent crystals can be formed either by closed packing of localized emergent particles, or by interacting macroscopic density waves with different wave vectors.
The proof of wave nature of emergent crytals relies on strict longrange order of the crystalline states. In real materials, existence of fluctuation, boundaries, defects, and microstructures all affects this precondition, for which an emergent particle in an emergent crystal may present particle properties locally. For instance, the L-shaped elongated skyrmions observed in supercooled Co 8 Zn 8 Mn 4 thin plates as a metastable phase 36 shows simultaneously a tendency towards the IPSQ phase, as well as the presence of local topological protection of the composing skyrmions. To deal with these effects, the present theoretical method needs to be developed. When the effect of thermal fluctuation is small enough to be regarded as a perturbation of the equilibrium state, the saddle-point method provides the first order approximation to account for the effect of fluctuation on the equilibrium properties 1,37 . Yet, when the long-range order is severely damaged, such as in the skyrmion glass state 38 , it will be extremely hard to construct an effective method based on modification of the present theory. Similar discussion can be made on the effect of boundary and sample size. We know that as the diameter of a nanodisk of chiral magnets increases, the magnetic state evolves from an isolated skyrmion to skyrmion clusters and then to SkX 39 . This example shows the competition between boundary effects and bulk effects. In magnetic states where the boundary effects dominate (i.e., isolated skyrmion and skyrmion clusters), the magnetic state depends sensitively on the shape and conditions of the boundary, for which the wave description may not be a good starting point towards an effective theory.

Discussion
We prove based on variational analysis that emergent crystals appearing in chiral magnets with long-range order has a wave nature, such that their exact mathematical expression is a Fourier series. Then we construct a systematic approach upon this analytical Fourier framework to obtain the equilibrium properties, the metastabiliy region in the phase diagram, and the phase transition path of any types of emergent crystals. We apply the method to analyze four types of emergent crystals, including the SkX, DiskX, the IPSQ phase, and the ferromagnetic phase. For SkX, we successfully explain its phase transition behaviors, such that it transforms to the IPSQ phase at low magnetic field and to the ferromagnetic phase at high magnetic field. We show that emergent crystals can be formed by interaction between spindensity waves with different wave vectors. This new formation mechanism differs emergent crystals from ordinary atomic crystals.

Methods
Free energy functional and rescaling process. We use the following free energy density functional to study magnetic skyrmions in cubic helimagnets where M denotes the magnetization vector, B denotes the external magnetic field, and T denotes the temperature. The five terms on the rhs. of Eq. (7) denote respectively the exchange energy density, the Dzyaloshinskii-Moriya (DM) interaction, the Zeeman energy density, and the second and fourth order Landau expansion terms. Eq. (7) can be reduced towðmÞ ¼ β K 2 wðmÞ in Eq. (1) by rescaling the spatial variables as To derive Eqs. (3)-(5), Eq. (1) should be described in the spherical coordinates asw Solving the equilibrium properties of an emergent crystal. Substituting the nth order Fourier representation of emergent crystals into Eq. (1), after integration we have where where q ij ¼ s i q ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi q 2 ijx þ q 2 ijy q and i ¼ ffiffiffiffiffiffi À1 p . The eigenvalues of A ij are (s i q − 1) 2 , 1 þ s 2 i q 2 , and (s i q + 1) 2 , and the corresponding unit eigenvectors can be solved as P ij1 ¼ 1 ffiffi 2 p s i q Àiq ijy ; iq ijx ; s i q h i T , P ij2 ¼ 1 s i q ½q ijx ; q ijy ; 0 T , P ij3 ¼ 1 ffiffi 2 p s i q ½iq ijy ; Àiq ijx ; s i q T . Hence m qij can be expanded as a linear combination of the three unit eigenvectors where c ijk ¼ c re ijk þ ic im ijk and c re ijk and c im ijk are real variables to be determined. Using Eq. (12), w per can be recasted as P n i¼1 P n i j¼1 c ij1 , which is positive semidefinite. By solving w per ¼ 0, we obtain a unique solution q = 1, |c ij2 | = |c ij3 | = 0 for all i, j, and |c ij1 | = 0 for i ≠ 1. This provides us with a valuable start-point to search for any possible emergent crystalline state in chiral magnets described by the free energy density functional given in Eq. (1). Taking 1st order Fourier representation for triangle lattice as an example, there are 22 independent variables to describe an emergent crystal, which can be written in a vector form as in which case the characteristic vector V C is replaced by the dominant vector V D (i.e., V D1 = v d U D + v 1 V R1 + v 2 V R2 + ... + v s V Rs , where V R1 , V R2 , ..., V Rs are unit vectors extracted from R 1 , R 2 , ..., R s ), and then do all the calculation again. By doing so, one derives the phase transition evolution path for an intrinsically unstable emergent crystal described by V C0 , and then the destination of this evolution: another emergent crystalline state that is metastable. In practice, one will find that there always exist two eigenvalues of D very close to zero for any solution V 0 which remain almost invariant with b. These two modes correspond to in-plane rigid body translation of the emergent crystals and should be distinguished from the real soft modes.

Data availability
The data that support the plots within this paper and other findings of this study are available from the author on reasonable request.