Skyrmion crystal formation and magnetic eld

The roles of magnetic ﬁeld and temperature in thermodynamic formation of skyrmion crystals (SkXs) have not been well-revealed to date. Here we present a uniﬁed theory about SkX formation and its fascinating thermodynamic behaviours. A chiral ﬁlm can have many metastable states with an arbitrary number of skyrmions up to a maximal value. A perpendicular magnetic ﬁeld makes a ﬁlm with Q m skyrmions the lowest energy state. Q m ﬁrst increases with the magnetic ﬁeld up to an optimal value and then decreases with the ﬁeld. The ﬁlm with the largest Q m at the optimal ﬁeld is an SkX. Outside of a ﬁeld window, states consisting of various stripes with low skyrmion number densities are thermal equilibrium phases while an SkX is metastable. Within the ﬁeld window, SkXs are the thermal equilibrium states below the Curie temperature. However, the time to reach an SkX state from a stripy phase would be too long at a low temperature. This causes a widely spread false belief that SkXs are metastable and stripy states are thermal equilibrium phase at low temperature and at the optimal ﬁeld. Our theory opens a new avenue for SkX manipulation and skyrmion-based applications.

SkXs form thermodynamically only under the assistance of an optimal magnetic field and near the Curie temperature. Once formed, however, SkXs can be metastable in very large temperature-field regions. For example, an SkX can exist at zero magnetic field and a field much higher than the optimal value by cooling an SkX first at the optimal field to low temperature followed by removal or rise of the magnetic field. An SkX disappears during the zero field warming and high field warming [11,18,19]. An SkX is believed to be a thermal equilibrium state only in a very narrow magnetic-field range and near the Curie temperature. Outside of the range, SkXs can only be metastable states that are far from equilibrium. At zero field and the temperature much lower than the Curie temperature, the equilibrium phase is believed to be a collection of stripe spin textures known as a helical state [20,21].
The thermal equilibrium state of a system comes from the competition between the internal energy and entropy. Entropy dominates the higher temperature phase while internal energy determines lower-temperature one. The entropy of a less-organized structure such as helical state is larger than a well-organized structure such as an SkX. Many of current views about SkXs are not consistent with this principles [11]. For example, as a thermal equilibrium state near the Curie temperature and at an optimal field, an SkX of low entropy structure must have a lower energy than a helical phase of relative high entropy structure. This putative conclusion contradicts the popular belief that the SkX becomes a metastable at a lower temperature and at the optimal magnetic field while a helical phase suddenly takes over the position of the thermal equilibrium state under the same conditions from the SkX. Obviously, the strange thermodynamic-pathdependent behaviours of SkXs and stripy phases, as well as SkX formation have not been satisfactorily explained yet.
In this letter, we show that a chiral magnetic film with an arbitrary number of skyrmions up to a critical value is metastable when skyrmion formation energy is negative and skyrmions are stripes. The energy and morphology of theses metastable states depend on the skyrmion number and the magnetic field perpendicular to the film. At zero field, the energy increases with skyrmion number or skyrmion number density. Thus thermal equilibrium phase below the Curie temperature at zero field should be helical states consisting of a few stripe skyrmions. At non-zero field, the film with Q m skyrmions has the lowest energy. Q m first increases with the field up to an optimal value then decreases with the field. Q m near the optimal field is large enough such that the average distance between two nearby skyrmions is comparable with skyrmion stripe width, and skyrmions form an SkX. Although the thermal equilibrium state should be an SkX near the optimal field and below the Curie temperature, an SkX appears thermodynamically only near the Curie temperature, and it takes too long for the system to reach this equilibrium state from a helical state at low temperature due to topological protection.
We consider here a magnetic thin film of thickness d in xy plane. The magnetic energy of the film is where A, D, K, H, M s , and µ 0 are ferromagnetic exchange stiffness constant, DMI coefficient, perpendicular easy-axis anisotropy, perpendicular magnetic field, the saturation magnetization, and the vacuum permeability, respectively; and m is magnetization unit vector. E = 0 is assumed for ferromagnetic state of m z = 1. It is known that isolated circular skyrmions are metastable state of energy 8πAd √ 1 − κ when κ = π 2 D 2 /(16AK) < 1 [22].

Results
Metastable structures under magnetic fields.  Figure 1(d1) shows how total Q (the left y−axis and the red curves) and energy E (the right y−axis and the blue curves) change with time. Within a very short time of order of sub-picoseconds, Q reaches its final stable values while it takes nanoseconds for the energy to decrease to their minimal values. The negative skyrmion formation energy explains well the stripe morphology of skyrmions that try to fill up the whole film in order to lower its energy, in contrast to circular skyrmions for positive energy [22]. Unexpectedly, the film can host an arbitrary number of skyrmions up to a large µ 0 H=0.4T µ 0 H=0.1T value (see also Fig. 2 below). At a low skyrmion number of Q = 10, the film is in a helical state with ramified stripe skyrmions of well-defined width of 8.1nm [23]. The film is in a helical state with rectangular stripe skyrmion at Q = 50 while it is an SkX of triangular lattice at Q = 150. When Q = 150, or skyrmion number density, is so large such that the distance between two neighbouring skyrmions is comparable to the stripe width, skyrmionskyrmion repulsion compress each skyrmions into a disklike object. Skyrmions at high skyrmion number density prefer a triangular lattice as shown in (a3) instead of the initial square lattice. Figure 1(b1-b3 and c1-c3) plot the metastable structures under fields µ 0 H = 0.1T (b1-b3) and µ 0 H = 0.4T (c1-c3) with the same initial states as those in (a1-a3). Stripes of m z > 0 (the grey regions) parallel to H expand while those of m z < 0 (the white regions) anti-parallel to H shrink. This is showed in Fig. 1 (b1, b2, c1, and c2) and spin profiles in the insets (see Supplementary Materials). Moreover, the amount of increase and decrease of white and grey stripes are not symmetric such that skyrmion-skyrmion repulsion is enhanced by the field, and SkXs tend to occur at lower skyrmion number density. This trend can be clearly seen in (c2) and (c3). Fig.  1(d2) and (d3) show similar behaviour for µ 0 H = 0.1T (d2) and 0.4T (d3) as their counterparts (d1) at zero field: Skyrmion number Q grows to final values rapidly in sub-picoseconds and E monotonically decreases to their minimum in nanoseconds. Figure 1 shows unexpectedly that the nature shape of skyrmions are various types of stripes when κ > 1 and at low skyrmion number density, in contrast to the current belief of all skyrmions are circular. We found also that one uncleation domain develops into one skyrmion.
Skyrmion-number dependence of energy. We compute below the skyrmion-number dependence of energy at a given magnetic field in order to understand the role of a field in SkX formation. We limit our study here to the Q far below its maximal value of more than 500. Physics around the maximal Q is not our current concern, and it will be investigated in future. Q nucleation domains of 5nm in diameter each arranged into a periodical lattice is used as the initial configuration to generate a steady structure of Q skyrmions with energy E. The left panel of Fig. 2 is the numerical results for µ 0 H = 0, 0.1, 0.2, 0.3, 0.39, 0.5, 0.6 and 0.7T. At zero field, E increases monotonically with skyrmion number. Thus, states with few skyrmions or low skyrmion number density are preferred. Since long stripe skyrmions have more way to deform than disk-like skyrmions, the entropy of a helical phase is larger than an SkX. Thus, helical states should always be the thermal equilibrium phase below the Curie temperature, and an SkX can be a metastable state at most. E is not sensitive to Q, thus the thermal equilibrium helical states can have various forms or morphology and contains different number of irregular skyrmions. Things are different when a magnetic field is applied. E of fixed Q increases with H, and E is minimal at Q m for a fixed field below a critical value. Q m first increases with H up to an optimal field of around µ 0 H = 0.3T and then decreases with H. Above µ 0 H = 0.7T (the stars and the brown), positive E means that ferromagnetic state of m z = 1 is the ground state with E = 0 and should be thermal equilibrium state below the Curie temperature when µ 0 H > 0.7T. Q m , at the optimal field of µ 0 H = 0.3T, can be as large as more than 141 or a skyrmion number density more than 3500/µm 2 at which two nearby skyrmions are in contact.
All skyrmions are compressed into a disk-like object and form an SkX. Strictly speaking, they are not circular, as evident from our simulations. All our results agree quantitatively with SkX experiments on MnSi [2,3].
Thermal effects and thermodynamic-path dependent phases. To substantiate our assertion that topological protection prevent a metastable helical state from transforming into a thermal equilibrium SkX state at µ 0 H = 0.3T and far below the Curie temperature, we study stochastic LLG at a finite temperature starting from the structure of Fig. 1(c1) with 10 skyrmions. We run the MuMax3 [24] at 20K and 29K, respectively below and near the Curie temperature of T C = 33K. Fig. 3(a) and (b) show the structures after 30ns evolution. Thermal fluctuations at 20K are not strong enough to neither create enough nucleation centres nor cut a short stripe into two that are the process of breaking the conservation of skyrmion number. The system is still in a helical state with Q = 11 skyrmions, one more than the initial value (see Movie 1 of Supplementary Materials). However, at 29K, a helical state transforms to the final thermal equilibrium SkX state with Q = 132 skyrmions through cutting stripes into smaller pieces and creating nucleation centers to generate more skyrmions thermodynamically (see Movie 2 of Supplementary Materials). For a comparison, we have also started from the SkX shown in Fig.  1(c3) Fig. 1(c1). (c) Structures after 10ns evolution at µ0H = 0.3T and T = 29K, starting from the initial configuration of Fig. 1(c3). (d,e) Structures after 30ns zero field warming (d) and warming (e), starting from the SkX in Fig. 1(c3). (f) An SkX of 133 skyrmions at µ0H = 0.3T and 29K comes from the thermally generated nucleation centres.
as shown in Fig. 3(c) after 10ns evolution (see Movie 3 of Supplementary Materials), 14 less than the starting value. This is expected because the average skyrmion number at thermal equilibrium should be smaller than Q m = 141 according to the energy curve of µ 0 H = 0.3T in Fig. 2. We demonstrate below that an SkX at zero field is not a thermal equilibrium state by showing the disappearance of an SkX in both zero field cooling and warming. Starting from the SkX in Fig. 1(c3) and gradually increasing (decreasing) the temperature from 0K (30K) to 30K (0K) at sweep rate of 1K/ns with H = 0 (see Movies 4 and 5 of Supplementary Materials), final structures shown in Fig.  3(d) (zero field cooling) and (e) (zero field warming) are helical state consisting of stripe skyrmions. In contrast, field cooling at the optimal field of µ 0 H = 0.3T from 30K to 0K at same sweep rate does not change the nature of the SkX. This is consistent with our assertion that SkXs are the thermal equilibrium states at µ 0 H = 0.3T below the Curie temperature.
The nucleation centres can also be thermally generated near the Curie temperature such that skyrmions can develop from these thermally generated nucleation centres, rather than from artificially created nucleation domains. To substantiate this claim, we carried out a MuMax3 simulation at 29K under the perpendicular field of µ 0 H = 0.3T, Fig. 3(f) is a snapshot of spin structure of the thermal equilibrium state for the same film size and with the same model parameters as those in other figures. An SkX with 133 skyrmions is observed. The birth of these skyrmions and how the SkX is formed can be seen from Movie 6 of Supplementary Materials. The tiny energy gain or loss from the transition between two states of different Q makes such a transformation difficult because of the conservation of skyrmion number under continuous spin structure deformation. This study shows that, similar to liquid drop formation, new skyrmions can be generated only from nucleation centres or by splitting a stripe skyrmion into two. These process require external energy sources such as the thermal bath and result in topological protection and energy barrier between states of different skyrmion numbers. Although the energy of Q = 141 SkX has the lowest energy at µ 0 H = 0.3T at zero temperature (Fig. 2), an initial state with a few stripe skyrmions would not assume the lowest energy state at low temperature ( Fig. 3(a)). This demonstrates the multiple metastable states of various Q and topological protection to prevent SkXs and helical states from relaxing to the thermal equilibrium phases. This new understanding can perfectly explain those fascinating appearance and disappearance of SkX and helical states along different thermodynamic paths [18,19]. For example, the persistence of an SkX stability in the field cooling at the optimal field is because the SkX is the thermal equilibrium state below the Curie temperature, not a metastable state far from the equilibrium as common wisdom believes. The disappearance of an SkX in zero field warming and high field warming is because helical state is the thermal equilibrium phase. At high enough temperature below the Curie temperature, thermal fluctuations can spontaneously generate enough nucleation centres such that the system can change its skyrmion number and reach its thermal equilibrium phase of either helical state of low skyrmion number density or SkX state of high skyrmion number density.

Discussion
Various stripes, including curved, ramified, maze, have long been observed in experiments and simulations [27][28][29][30][31][32][33][34], but there is no good description about those complex spin structures to date. For short race-track-like stripes, a notion of meron with ½ skyrmion number or bimeron [34][35][36] were used to describe one end or whole structure. This local description is not accurate and not necessary, an those bimerons should be correctly called skyrmions in order to reflect a holistic view of the spin structure. Furthermore, the claim of each stripe end carrying ½ skyrmion number [35] is not correct. This claim cannot explain ramified stripe skyrmion with multiply ends. Topological charge distribution in irregular stripe skyrmioins is not what described in meron notion. In literature [20,21] spiral, helical, and conical states are used to describe various spin textures whose spins vary along one particular direction within a 2D or 3D materials. In our holistic description, many of them, if not all, are stripe skyrmions.
We have presented results for the interfcial DMI so far, similar results are also true for bulk DMI. The only difference is the change from Neel-type of stripe walls to the Bloch-type stripe walls [37]. Our findings, after firmly confirmed, can push the study of topological Hall effect to a new level because topological Hall effect can be investigated against much wide range of condensed skyrmion states.
Both helical states and SkXs are the collections of skyrmions. Stripes, not disks, are the natural shapes of skyrmions when their formation energy is negative. The distinct morphologies of helical and SkX states come from the skyrmion number density. Skyrmions become disk-like objects and in a triangles lattice due to the compression and strong repulsion among the highly packed skyrmions. Unexpectedly, there are enormous number of metastable states with an arbitrary number of skyrmions up to a critical value of above 500 at zero temperature. The physics around this value needs further studies. The energy of these states depends on the skyrmion number and the magnetic field. The role of a magnetic field in SkX formation is to create the lowest energy state being a skyrmion condensate with a high number of skyrmions. Our findings have profound implications in skyrmionbased applications.

Methods
Spin dynamics in a magnetic field is governed by the Landau-Lifshitz-Gilbert (LLG) equation, where γ and α are respectively gyromagnetic ratio and Gilbert damping constant. H eff = 2A µ0Ms ∇ 2 m+ 2K µ0Ms m zẑ + Hẑ + H d + H DM + h is the effective field including the exchange field, the anisotropy field, the external magnetic field alongẑ, the demagnetizing field H d , the DMI field H DM , and a temperature-induced random magnetic field of magnitude h = 2αk B T /(M s µ 0 γ∆V ∆t), where ∆V , ∆t, and T are the cell volume, time step, and the temperature, respectively [24,38].
In the absence of energy sources such as an electric current and the heat bath, the LLG equation describes a dissipative system whose energy can only decrease [25,26]. Thus, solving the LLG equation is an efficient way to find the stable spin textures of Eq. (1). Note that finding metastable states of a given physical system does not need to include thermal effect. In this study, we choose A = 0.4pJ/m, D = −0.33mJ/m 2 , K = 0.02MJ/m 3 , M s = 0.15MA/m to mimic chiral magnetic films such as MnSi [39,40] that satisfies κ > 1 and supports stripe skyrmions [37]. Periodic boundary conditions are used to eliminate boundary effects and the MuMax3 package [24] is employed to numerically solve the LLG equation with mesh size of 1 nm× 1 nm× 1 nm. It should be mentioned that the physics discussed here does not depend on the boundary conditions. However, different boundary conditions have different confinement potentials (or potential well depth) for skyrmions, and can affect the maximal skyrmion number density below which a condensed skyrmion state is metastable. The number of stable states and their structures should not depend on the Gilbert damping constant. We use a large α = 0.25 to speed up our simulations.

Data availability.
The data that support the plots within this paper and other findings of this study are available from the corresponding author on reasonable request.
References H.T.W. contribute equally. All authors discussed the results and commented on the manuscript.

Additional information
Supplementary Information accompanies this paper at http:// Competing interests: The authors declare no competing interests. X. R. Wang is an Editorial Board Member for Communications Physics, but was not involved in the editorial review of, or the decision to publish this article.