Stripe skyrmions and skyrmion crystals

Skyrmions are important in topological quantum field theory for being soliton solutions of a nonlinear sigma model and in information technology for their attractive applications. Skyrmions are believed to be circular and stripy spin textures appeared in the vicinity of skyrmion crystals are termed spiral, helical, and cycloid spin orders, but not skyrmions. Here we present convincing evidences showing that those stripy spin textures are skyrmions,"siblings"of circular skyrmions in skyrmion crystals and"cousins"of isolated circular skyrmions. Specifically, isolated skyrmions are excitations when skyrmion formation energy is positive. The skyrmion morphologies are various stripy structures when the ground states of chiral magnetic films are skyrmions. The density of skyrmion number determines the morphology of condensed skyrmion states. At the extreme of one skyrmion in the whole sample, the skyrmion is a ramified stripe. As the skyrmion number density increases, individual skyrmion shapes gradually change from ramified stripes to rectangular stripes, and eventually to disk-like objects. At a low skyrmion number density, the natural width of stripes is proportional to the ratio between the exchange stiffness constant and Dzyaloshinskii-Moriya interaction coefficient. At a high skyrmion number density, skyrmion crystals are the preferred states. Our findings reveal the nature and properties of stripy spin texture, and open a new avenue for manipulating skyrmions, especially condensed skyrmions such as skyrmion crystals.

Skyrmions, originally used to describe resonance states of baryons [1], were first unambiguously observed in the form of skyrmion crystals in various chiral magnets and by various experimental techniques [2][3][4][5][6][7]. Magnetic skyrmions are topological non-trivial spin textures of magnetic films characterized by a skyrmion number of Q = 1 4π m · (∂ x m × ∂ y m)dxdy, here m is the unit vector of magnetization. Q must be an integer for an infinite magnetic film, and a non-zero Q spin texture is called a skyrmion of skyrmion number Q. It was isolated circular skyrmions, not a skyrmion crystal, that were predicted in early theories [8,9]. Isolated skyrmions were indeed observed later in confined structures and films [10][11][12][13][14][15]. It is an experimental fact that skyrmion crystals form in very narrow magnetic-field-temperature windows. Outside of the windows, stripy phases appear. The stripy phases, which can even coexist with skyrmion crystals, are in fact easier to form than a skyrmion crystal does. In contrast, those stripy phases do not appear together with isolated skyrmions. Interestingly, the stripy phase was observed many years before the observations of skyrmion crystals [16]. These stripy phases are called helical, spiral, and cycloid spin orders. A one-dimensional model [9,[17][18][19][20][21] was used to describe the rotation of spins perpendicular to stripes. To date, a holistic description of various stripy structures, especially ramified stripes and stripy maze, [22][23][24][25][26][27][28][29] is lacking. The general belief is that those stripes are not skyrmions and have zero skyrmion numbers [30] although some race-track stripes are called merons with 1/2 skyrmion number or bimerons [27,31,32].
In this letter we show that stripy magnetic textures appear in a chiral magnetic film with Dzyaloshinskii-Moriya interaction (DMI) are actually irregular skyrmions. Each stripy texture has exactly skyrmion number 1. The irregular shape is due to the negative skyrmion formation energy (relative to ferromagnetic state) when the ferromagnetic state is not the ground state. For a given system, the morphology of these skyrmions are random when the skyrmion number density is low. At extremely low density, magnetic textures are ramified stripes. The exact appearance of each pattern is very sensitive to the initial spin configuration and its dynamical path. The basic building blocks of irregular random skyrmions are stripes of well-defined width. The optimal width comes from the competition between the Heisenberg's exchange energy and the DMI energy that respectively prefer a larger and smaller width in order to minimize exchange energy cost and maximize the negative skyrmion formation energy gain. Unexpectedly, this exchange energy and DMI energy dependence of width is opposite to the skyrmion size of an isolated skyrmion that increases with DMI interaction and decreases with exchange energy [34].
We consider an ultra-thin ferromagnetic film of thickness d in the xy plane. The film has an exchange energy E ex with exchange stiffness constant A, an interfacial DMI energy E DM with DMI coefficient D, an anisotropy energy E an with a perpendicular easy-axis anisotropy K, and the Zeeman energy E Ze in a perpendicular magnetic field H. The total energy E total reads where M s is the saturation magnetization and µ 0 is the vacuum permeability. The integration is over the whole film. The energy is set to zero, E total = 0, for the ferromagnetic state of m z = 1. We have assumed m is uniform in thickness direction. The demagnetization effect is included in E an through the effective anisotropy K = K u − µ 0 M 2 s /2 corrected by the shape anisotropy, here K u is the perpendicular magnetocrystalline anisotropy. This is a good approximation when the film thickness d is much smaller than the exchange length [34]. It is known that isolated circular skyrmions are metastable state of energy 8πAd √ 1 − κ when κ = π 2 D 2 /(16AK) < 1 [34]. Here we use MuMax3 simulator [35] to numerically solve the Landau-Lifshitz-Gilbert (LLG) equation for the stable states in the opposite regime of κ > 1 where circular skyrmions are not stable states [34], see Method.
For a sample of 300nm × 300nm × 0.4nm under the periodical boundary conditions and with A = 10pJ/m, D = 6mJ/m 2 , K = 0.49MJ/m 3 , and M s = 0.58MA/m that are typical values for those chiral magnets supporting skyrmion crystals [30,36,37], κ > 1, so that single domain and isolated circular skyrmion are not stable any more [34]. We will start from a small nucleation magnetic domain with sharp domain wall to speed up skyrmion formation dynamics although stripe skyrmions will appear spontaneously due to thermal or other fluctuations in reality. Figure 1 shows how various initial states in the top panel (a1, b1, c1) evolve according to the LLG equation with α = 0.25. α will not change physics described below because we are interested in the spin textures of the LLG equation that do not vary with time. However, α can change the evolution path and energy dissipation rate so that it can change the intermediate states shown in (a2-c2) and influence which fixed point to pick when a system has many fixed points (textures) like the current situation. This is similar to the sensitiveness of attractor basins to the damping in a macro spin system [38]. The initial configurations are obtained by reversing spins in the white regions in (a1-c1) from m z = −1 to m z = 1 such that the configurations have very sharp domain walls and zero skyrmion number Q = 0. After a short time of an order of picosecond, the initial states transform into irregular structures of skyrmion number Q = 1, no matter whether the initial shape is circular or non-circular as shown in (a4-c4) where the time is in logarithmic scale. As time goes on, Q (the blue lines in a4-c4) stays at 1, and system energy E total (the red curves) is negative and keeps decreasing until it reaches a stable ramified stripy spin texture. Clearly, this irregular ramified spin texture is a non-circular skyrmion whose formation energy is negative. The negative formation energy explains why the 58MA/m for a sample of 300nm × 300nm × 0.4nm. (a1,b1,c1) are different initial configurations of a disk of diameter 20nm (a1), a hexagon of side length 10nm (b1) and a square of length 20nm (c1). (a2,b2,c2) are intermediate states at 0.3ns with irregular shapes due to the negative formation energy. (a3,b3,c3) are the final stable pattern with irregular ramified stripes. Skyrmion charge density ρ is encoded by colours (the blue for positive and the red for negative) while the gray-scale encodes mz. The dark black lines denote mz = 0. The positive and negative charges exist respectively only around convex and concave areas. The spin profile across the stripes at the green n ○ is considered. (a4, b4, c4) are the evolution of total energy E total and topological skyrmion number Q (t is in the logarithmic scale). Q reaches the skyrmion number 1 within 1ps. Clearly, the skyrmion number is a constant and the total energy is negative and approaches a constant almost independent from the initial configurations.
skyrmion prefers to stretch out to occupy the whole space to lower its energy. This process is clearly demonstrated in the movie in the Supplementary Material, showing how the system evolves from (a1) to (a3) and how Q grows from 0 to 1. The simulations show also that the exact pattern is very sensitive to the initial configuration and dynamical/kinetic paths as well as how energy dissipates. In a real system, the process shall proceed spontaneously from any randomly generated nucleation centre no matter by a thermal fluctuation or an intentional agitation.
We use color to encode the skyrmion charge density ρ defined as ρ = m · (∂ x m × ∂ y m)/(4π) as indicated by the color bar in Fig. 1. Interestingly, ρ is non-zero only around convex (positive) and concave (negative) areas. ρ is almost zero along the long straight stripe. This may explain why the skyrmion number of stripy textures were thought to be zero in the literature [30]. As shown in Fig.  1(a4-c4), the sum of all positive and negative skyrmion charge is always quantized to 1 protected by the topology. The continuous decrease of energy E total also indicate the final morphology of the ramified stripe skyrmions is not unique and depends on dynamical path of the system evolution that in turn relates to the initial configuration.
One striking feature about the stripes shown in Fig.  1 is a well-defined stripe width that is usually referred to a given wave vector in experiments. The competition between the exchange energy and DMI energy can easily lead to a A/D dependence [9,[17][18][19][20]. However, it is highly non-trivial to understand the effect of the magnetic anisotropy on stripe width. In order to understand the underlying physics of this well-defined width, one needs a good spin texture profile along the direction perpendicular (x) to the stripes. We found that black (m z ≤ 0) and white (m z ≥ 0) stripes can be approximated by Θ(x) = 2 arctan sinh(L/2w) sinh(|x|/w) and Θ(x) = 2 arctan sinh(|x|/w) sinh(L/2w) (|x| ≤ L/2), respectively. Θ is the polar angle of the magnetization at position x and x = 0 is the centre of a stripe. L and w measure respectively the stripe width and skyrmion wall thickness as schematically illustrated in Fig. 2(a). Figure 2(b) demonstrates the excellence of this approximate profile for a set of model parameters of A = 10pJ/m, D = 6mJ/m 2 , K = 0.49MJ/m 3 , and M s = 0.58MA/m. The y−axis is m z and x = 0 is the stripe centre where m z = 1. Different symbols are numerical data across different stripes labelled by the green n ○ in Fig. 1(a3-c3). The solid curve is the fit of profile cos Θ(x) with L = 10.78nm and w = 3.05nm. All data from different stripes falling onto the same curve demonstrates that stripes, building blocks of pattern, are identical.
Using the excellent spin profile, one can obtain magnetic energy of a film filled with the stripe skyrmions as a function of L and w. Minimizing the energy against L and w allows us to obtain A, D, and K dependence of stripe width L and skyrmion wall thickness w. In terms of L, = L/(2w), ξ = A/D, and κ, the total system energy density of a film filled by such stripes is where [sinh 2 ( )+cosh 2 ( x)] 2 dx, and g 2 ( ) = 1 0 [sinh 2 ( )+cosh 2 ( x)] 2 dx. The optimal stripe width obtained from minimizing energy E is L = a 4A πD , where a depends weakly on κ = π 2 D 2 /(16AK) > 1. The physics of this result is clear: DMI energy is negative, and one can add more stripes by reducing L such that the total energy will be lowered. On the other hand, the exchange energy will increase with the decrease of L. As a result, L is proportional to A/D that is opposite to the behaviour of size of an isolated skyrmion whose size increases with D and decreases with A [30,34]. These theoretical results agree very well with micro-magnetic simulations as shown in Fig. 2(c) with a = 7.61 for κ = 1.3; a = 6.82 for κ = 1.8; a = 6.47 for κ = 2.9. The dependence of L on κ for a fixed A/D = 1nm is shown in the inset. One can see that L depends weakly on large κ 1. Unexpectedly, L are the same no matter whether we have one, two, or more ramified and non-ramified stripe skyrmions. This is reflected in the skyrmion number density independence of L as shown in Fig. 2(d).
To understand why condensed stripe skyrmions and skyrmion crystals can appear together in a given chiral magnetic film when the material parameters are fixed, we try to increase the number of skyrmions in our film of 300nm × 300nm × 0.4nm. Encouraged by the results in Fig. 1 that each nucleation domain creates one irregular stripe skyrmion, we place 2, 100 and 169 small disk domains of diameter 10nm as shown in Fig. 3(a1-c1) and let them to evolve according to LLG equation. Fig. 3(a2-c2) are the final steady states. As expected, we indeed obtained two irregular ramified stripe skyrmions (a2). For the case of 100 skyrmions, some of them have rectangular shape and are arranged in a nematic phase while the rest of skyrmions look like disks and are in a lattice structure (b2). In the case of 169 skyrmions, skyrmions are disklike and are in a triangular lattice (c2). The skyrmion nature of the spin textures can be clearly confirmed by the change of Q (blue) with time as shown in Fig. 3(d) (the solid, dashed, and dotted lines for Q = 2, 100, and 169 respectively). (d) shows also how system energy (the red lines) changes with time. One interesting feature is that total energy is not sensitive to skyrmion number density before skyrmion-skyrmion distance is comparable to the optimal stripe width and skyrmions take stripe shape. The system energy starts to increase with the skyrmion number density, and condensed skyrmions transform from rectangular stripe skyrmions in nematic phase into circular skyrmion crystal. This feature does not favour skyrmion crystal formation, and may explain why skyrmion crystals were observed in the presence of an external magnetic field: With only one skyrmion in the whole system as those shown in Fig. 1, the net magnetic  Fig. 1(a3-c3), are mz = cos Θ from MuMax3. The green solid line is the fit of Θ(x) = 2 arctan sinh(|x|/w) sinh(L/2w) . (Θ ≤ π/2) and Θ(x) = 2 arctan sinh(L/2w) sinh((|x−L|))/w) (Θ ≥ π/2) with L = 10.78nm and w = 3.05nm. x = 0 is the centres of white stripes. All data from different stripes fall onto the same curve means that stripes of the same width are basic building blocks of stripe skyrmions. (c) A/D-dependence of stripe width L for various κ = π 2 D 2 /(16AK) = 1.3 (the green stars); 1.8 (the blue squares); and 2.9 (the red circles). The solid lines are the fits of L = aA/D with a = 7.61 for κ = 1.3; a = 6.82 for κ = 1.8; a = 6.47 for κ = 2.9. Inset: the dependence of L on κ for A/D = 1nm. Symbols are the numerical data and the dashed line is the theoretical prediction without any fitting parameter. (d) The dependences of stripe width L on skyrmion number density for A = 10pJ/m, D = 6mJ/m 2 , K = 0.49MJ/m 3 , and Ms = 0.58MA/m. The number of skyrmions varies from 1 to 113 in our sample of 300nm × 300nm × 0.4nm under the periodical boundary conditions. The dash line is L = 10.52nm. moment is zero because spin up (white area) and down spin (black area) are almost equal so that magnetic moment cancel each other. As skyrmion number increases, the net magnetic moment increases and becomes nonzero.
We have showed that stripe and ramified stripe skyrmions are essentially the same as the circular skyrmions in skyrmion crystals. The difference in skyrmion shapes at different skyrmion number density come from the skyrmion-skyrmion interaction. When the average distance between two nearby skyrmions is order of the stripe width, the skyrmion-skyrmion repulsion compress skyrmions into circular objects. This understanding permits a skyrmion crystal in the absence of an external magnetic field as long as one can use other means to add more skyrmions into a film such as a scanning tunnelling tip [6].
We have presented results for the interfacial DMI so far, similar results are also true for bulk DMI. As shown in Fig. 4, one ramified stripe skyrmion has very similar structures for interfacial (a) and bulk (b) DMIs when we start from the same initial configuration with the same interaction strength. The only difference is the change of Neel-type of stripe wall to the Bloch-type stripe wall as shown by the color-coding in the figures.
A perpendicular magnetic field can modify the morphology and width of a stripe without changing its skyrmion number. Stripe width increases (decreases) To speed up the system reaching its final stable states, we use a large Gilbert damping constant of α = 1. (a1,b1,c1) are the initial configurations and (a2,b2,c2) are the final stable pattern. skyrmion charge density ρ is encoded in colours (the red and the blue for positive and negative ρ respectively) while mz is encoded in the grayscale. The dark black lines denote mz = 0. The positive and negative charges exist respectively only around canvex and concave areas. (d) is the evolution of total energy E total (the left y-axis) and topological skyrmion number Q (the right y-axis with three ranges). Clearly, Q reaches a constant integer in a very short time and Q = N = 2, 100 and 169 respectively. The total energy is negative and approaches a constant, indicating a stable state.
with the field when out-of-plane spin component is antiparallel (parallel) to the field.
It requires certain amount energy to destroy a stable or metastable state by definition. Thus, all skyrmions discussed in this paper are stable against thermal noise as long as the thermal energy is smaller than their potential barriers. In the Supplemental Material, we provide a video to show that the state in Fig. 3(b2) is stable at 50K.
Skyrmions provide a fertile ground for studying funda-  mental physics. For example, the topological Hall effect is a phenomenon about how non-collinear magnetization in skyrmion crystals affect electron transport. Knowing stripy phases are also condensed irregular skyrmions, it expands surely arena of topological Hall effect. We can not only study how electron transport be affected by the skyrmion crystals, but also skyrmions in other condensed phases such as nematic phases, or how the elongation and orientation of stripe skyrmion affect Hall transport. With the new discovery of stripe skyrmions, it will also allow us to investigate the interplay of topology, shape, spin and charge. One can investigate how the topology, the local and global geometries affect spin excitations separately. The assertion that ramified stripes and other stripes appeared together with skyrmion crystals are irregular skyrmions are firmly confirmed by the nano-magnetic simulations. These stripes have a well-defined width that is from the competition between exchange interaction energy and DMI energy. The inverse of the width is the well-known wave-vector used to describe the spiral spin order in the literature [16]. In contrast to isolated skyrmions whose size increases with DMI constant and decreases with the exchange stiffness constant, the stripe width increases with the exchange stiffness constant and decreases with the DMI constant. Counter-intuitively, skyrmion crystals are highly compressible like a gas, not like an atomic crystal. This detail property need a careful further study. We believe that our findings should have profound implications in skyrmion-based applications.

Methods
Numerical simulations. Spin dynamics is governed by the Landau-Lifshitz-Gilbert (LLG) equation, where m, γ, α are respectively the unit vector of the magnetization, gyromagnetic ratio, and the Gilbert damping. H eff = 2A∇ 2 m + 2Km zẑ + H d + H DM is the effective field including the exchange field characterized by the exchange stiffness A, crystalline anisotropy field, demagnetizing field H d , and DMI field H DM .
In the absence of energy source like an electric current, LLG equation describe a dissipative system whose energy can only decrease [39,40]. Thus, solving LLG equation is an efficient way to find the stable spin textures. In the present case, we apply periodic boundary conditions to eliminate the edge effects. We use the Mumax3 package [35] to numerically solve the LLG equation with mesh size of 1nm × 1nm × 0.4nm for those skyrmions of L > 5nm. For skyrmions of L < 5nm , the mesh size is 0.1nm × 0.1nm × 0.4nm. The number of stable states and their structures should not depend on the Gilbert damping constant, but spin dynamics is very sensitive to α. To speed up our simulations, we use large α of 0.25 and 1. We consider only material parameters that supports condensed skyrmion states. The skyrmion size L is obtained directly from numerical data or by fitting the simulated spin profile to Θ(x) = 2 arctan sinh(L/2w) sinh(|x|/w) (the black stripes for m z ≤ 0) and Θ(x) = 2 arctan sinh(|x|/w) sinh(L/2w) (white stripes for m z ≥ 0) with stripe width L and wall width w. Here −L/2 ≤ x ≤ L/2 and x = 0 is the centre of a stripe. It should be pointed out 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 stripe skyrmions, and can affect the maximal skyrmion number density below which a condensed skyrmion state is metastable.
Energy density of a film filled by stripes.
Following similar assumptions as those in Ref. [19], we can computer the energy density of a system filled with stripes of width L that parallel to the y−axis and are periodically arranged along the x−axis. If Θ is the polar angle of spins and spin profile is Θ(x) = 2 arctan sinh((|x−nL|)/w) sinh(L/2w) + nπ with integer n and x ∈ ((n − 0.5)L, (n + 0.5)L) as shown in Fig. 2(a), the total energy of the film in the range of y 1 < y < y 2 and In terms of L, , ξ = A/D and κ = π 2 D 2 /(16AK), we obtain Eq. (2), E = 4D ξ ξ 2 L 2 g 1 ( ) − πξ 4L + π 2 64κ g 2 ( ) .
The first and the third terms on the right hand side are positive, thus To have negative E, κ must be larger than √ g 1 g 2 ≥ 1.
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.