New spiral state and skyrmion lattice in 3D model of chiral magnets

We present the phase diagram of magnetic states for films of isotropic chiral magnets calculated as function of applied magnetic field and thickness of the film. We have found a novel magnetic state driven by the natural confinement of the crystal, localized at the surface and stacked on top of the conical bulk phase. This magnetic surface state has a three-dimensional (3D) chiral spin-texture described by the superposition of helical and cycloidal spin spirals. This surface state exists for a large range of applied magnetic fields and for any film thickness beyond a critical one. We also identified the whole thickness and field range for which the skyrmion lattice becomes the ground state of the system. Below a certain critical thickness the surface state and bulk conical phase are suppressed in favor of the skyrmion lattice. Unraveling of those phases and the construction of the phase diagram became possible using advanced computational techniques for direct energy minimization applied to a basic 3D model for chiral magnets. Presented results provide a comprehensive theoretical description for those effects already observed in experiments on thin films of chiral magnets, predict new effects important for applications and open perspectives for experimental studies of such systems.


Introduction
Chiral magnets (ChM) are a distinct class of magnetic crystals. Contrary to the classical ferroand antiferromagnets the ground state of ChM is an incommensurate homochiral spin spiralthe spiral with an unique sense of magnetization rotation, see figure 1(a). The interaction which is responsible for the stabilization of such spin spirals is an antisymmetric exchange, also known as Dzyaloshinskii-Moriya interaction (DMI) [1,2], which appears in crystals with broken inversion symmetry. For instance, to this class of magnetic materials belong different Si-and Ge-based alloys such as MnSi [3] and FeGe [4], Mn 1−x Fe x Ge [5], Mn 1−x Fe x Si [6], Fe 1−x Co x Si [7]. Such alloys with B20 crystal symmetry can be referred to as the distinct class of so-called isotropic ChM (IChM). Such a classification reflects the dominant role of DMI and Heisenberg exchange, which are assumed to be isotropic in all spatial directions, while in the frame of the basic model, the contribution of magnetocrystalline anisotropy can be neglected.
The special interest to such materials arose after the breakthrough results on the direct observation of magnetic skyrmions in thin films of ChM [7]. The experimental discovery of magnetic skyrmions together with the conceptual idea of a revolutionary new type of magnetic memory gave an additional impetus to the research in this field [9,10].
Here we present results of our theoretical calculations for the magnetic field induced transitions in thin films of IChM, which allowed us to identify the critical film thickness above which the skyrmions never appear as the ground state of the system. Realizing that even in the simplest one-dimensional (1D) model of ChM the ground state exhibits periodic modulations [11], in the three-dimensional (3D) systems one should expect the appearance and coexistence of complex inhomogeneous phases, which are of major interest for both fundamental research and practical applications. Indeed, in this paper we identify an earlier unknown phase, that we term stacked spin spirals phase (StSS). The calculated phase diagram exhibits a wide range of existence of this new phase. This new phase exhibits periodic modulations in all three spatial directions and can be considered as the coexistence of the conical phase and the complex spin spirals localized near the surfaces. We found that such a state should appear in a wide range of film thicknesses as well as in the bulk samples. We provide a comprehensive description of this phase together with our suggestions how such a state can be experimentally detected.  The vector field and cross sections of corresponding isosurfaces, n z = 0, for: skyrmion tube with homogeneous magnetization along the radial symmetry axis (d), for skyrmion tube in thin film of chiral magnet with inhomogeneous magnetization in all three spatial directions and the twist induced by the free surface (e), chiral bobbers -hybrid particlelike states localized near the surface of the chiral magnet characterized by smooth magnetization distribution and presence of singularity -chiral Bloch point at finite distance from the surface (f), for details see Ref. [19]. Note, the section of the isosurfaces in (d)-(f) corresponds to n z > 0, n x < 0.
It is worth mentioning that the existence of the magnetic skyrmions in magnetic crystals with DMI has been theoretically predicted by Bogdanov in 1989 [8], but it took more then twenty years to discover them experimentally. One of the reason behind is the modification of the conventional energy balance in thin films of ChMs in comparison to the bulk crystals. Indeed, according to the theory developed by Bogdanov and co-workers, see e.g. Ref. [12], the chiral skyrmion tubes, see figure 1(d), in bulk crystals of ChMs with relatively weak magnetocrystalline anisotropy may appear only as a metastable state. Their energy is always higher than the energy of the conical phase, figure 1(b), which dominates in bulk crystals almost in the whole range of magnetic fields and temperatures. This result is consistent with many experimental studies on bulk ChM over the past decades, where skyrmions were not observed. The only exception is the so-called A-phase -the high-temperature region in the bulk phase diagram, just below the ordering temperature, where skyrmions has been proposed to exist due to thermal fluctuations [13]. However, the real nature of the A-phase still is under debate [13,14,15]. The theoretical models predict that in bulk crystals of ChMs, the magnetic skyrmions can be stabilized due to the strong cubic or uniaxial anisotropy [12,16] or special crystal symmetry, which suppress the formation of the conical phase [8,17]. However, such theoretical models as well as thermal fluctuations can not be considered as the main mechanism for skyrmion stabilization in stand-alone thin films of B20alloys where the skyrmions are observed in a wide range of temperatures, much lower than the ordering temperature.
As has been shown in Ref. [18] the key to the understanding of the mechanism of skyrmion stabilization in thin films of ChM are presence of the free surfaces and the three-dimensional rather than the two-dimensional structure of the equilibrium skyrmions. As shown in figure 1(e), the solution for the skyrmion tube in a thin film is characterized by the twist of magnetization with respect to the normal vector of the film surfaces. One can compare the twisted iso-surface of the skyrmion tube in figure 1(e) and the homogeneous non-twisted skyrmion tube in (d). The magnetic moments in the top surface layer are slightly turned towards the center of the skyrmion, while in the bottom of the film they are slightly turned outwards the center. The spin structure on the top and bottom surfaces corresponds to a certain intermediate configuration in between of pure Bloch-and Neel-type of skyrmions. Such surface induced twist propagates from one surface to another through the whole film. Note, it mainly affects the magnetic spins near the surfaces where spins are weakly coupled to each other because of reduced number of neighbors at the free surface. Far from the film surface, the spin structure of such a 3D skyrmion remains almost the same as in a homogeneous skyrmion tube.
The energy gain by the DMI contribution accumulated along the film thickness reduces the total energy of the state such that within a certain range of magnetic fields and film thicknesses the skyrmion tube becomes energetically more favorable than the conical phase. Moreover, it has been shown earlier that the same effect is responsible for the stability of chiral bobberparticlelike objects localized in all three spatial directions near the surface of the chiral magnet [19], see figure 1(f). Note, because of the presence of a singularity, the chiral Bloch point, at low temperatures the chiral bobbers appear exclusively as a metastable state.
It is obvious that for very thick films, the relative energy contribution of the surface twist becomes very small, while the main contribution to the energy of skyrmion comes from the volume part of the film. Thereby, there should be a critical thickness above which the energy gain of the surface twist is not anymore sufficient to provide enough energy gain to stabilize a skyrmion lattice. In order to identify the range of thicknesses and magnetic fields defining the range of stability for the skyrmion lattice and other states, we have calculated a phase diagram for the film of isotropic ChM in a wide range of thicknesses and applied magnetic fields, which is presented in section 3.1.

Model
The basic model for IChM include three main energy terms: the Heisenberg exchange interaction, the DMI and Zeeman energy term [23,24]: where n ≡ n(r) is a continuous unit vector field defined everywhere except at the singular points. E 0 is the energy of the saturated ferromagnetic state. A and D are micromagnetic constants for exchange and DMI, respectively, and M s is the magnetization of the material -the total magnetic dipole moment per unit volume.
We use the continuum model as the most general approach to describe long-period incommensurate magnetic structures. The results presented here can be easily generalized for a wide class of the systems. The functional (1) has to be considered as a continuum limit for the classical spin models e.g. as simplified models considering a simple cubic lattice [25], and advanced models, which take into account the exact B20 crystal symmetry [26,27] .
The lowest period of an incommensurate spin spiral and equilibrium period of the conical phase, L D , as well as the critical field corresponding to the saturation field of the conical phase, H D , have analytic solutions, which couple the material parameters with experimentally measurable quantities [11,17]: The comparison of the energy density of each of the equilibrium states obtained by direct energy minimization of the functional (1) allows one to identify the geometrical and material parameters corresponding to the phase transitions. For details of the energy minimization technique and calculation of the phase diagram, see section 4.

Phase diagram
In the phase diagram presented in figure 2(a) the thickness of the film, L, and magnetic field H are given in reduce units, where L D and H D are the functions of material parameters A, D and M s , see Eq. (2). The unique pair of parameters L D and H D can be considered as a fingerprint of each particular IChM. They can be measured experimentally, which allows to rescale the phase diagram, figure 2(a), in real units of film thickness, L and magnetic field H.
The solid lines in figure 2(a) correspond to the first-order phase transitions between helical spin spiral (red), skyrmion lattice (grey), conical phase (white) and new up-to-now unknown phase, which we call stacked spin-spirals state (yellow) and is discussed in detail in section 3.2. The horizontal dashed line indicates the second order phase transition between conical and saturated ferromagnetic state (blue).
The phase diagram for a thin layer of IChM for relatively narrow range of thicknesses, 0< L/L D ≤4, had been presented earlier in Ref. [19] . An attempt to generalize such a phase diagram for a wider range of thicknesses, 1≤L/L D <∞, had been announced recently [20]. The phase diagram presented in Ref. [20] contains four distinct states and one triple-point. We have tested these data and got qualitatively and quantitatively different results. The phase diagram presented in this contribution contains five distinct states and two triple points.
The reflects the relative contribution of the surface induced twists, which increases with decreasing thickness. There is another critical point, for L/L D < 0.68, the conical phase is totally suppressed and becomes energetically unfavorable in the whole range of fields.
It is important to mention that the chiral surface twist discussed above also introduces an additional modulation in the helical spiral state. Note, the k-vector of such helical spiral lies in the plane of the film, orthogonal to the applied field. The surface induced modulation reduces the energy of the helix and in a certain range of fields makes it energetically more favorable than the conical state. Such a behavior of the system is totally different to the one of the bulk crystals, where theoretically any infinitesimal magnetic field leads to convergence of the helix to the conical phase or more precisely to the StSS according to results presented here, see right panel in figure 2(a).
Finally, we have to emphasize that the effect of the chiral surface twist is not restricted to the film surfaces, but also appears on the side edges of the sample. The presence of the edge twist effect has been reported earlier in Ref. [21] and has been confirmed recently by directly observation with the Lorentz Transmission Electron Microscopy (TEM) [22].
It is worth to emphasize that, the continuum 3D model of IChM strictly speaking does not converge to a simple 2D model even in the limiting case of L/L D → 0. In order to illustrate such a discrepancy we added the left panel in figure 2(a), which corresponds to the phase diagram of 2D IChM valid for the single isotropic monolayer or multilayer with interface induced DMI and for ChM of particular crystal symmetry e.g. C nv , D 2d and S 4 [8,17]. A wide range of the phase diagram is occupied by the newly found stacked spin spirals state (StSS). The triple point II defines the limiting thickness above which the StSS may appear as the global energy minimum.  We found that the relative orientation of the propagation vectors of the spin spirals on the top and bottom surfaces, k t and k b , respectively, is thickness dependent. In other words, the angle β tb between k t and k b for the equilibrium StSS varies as function of the film thickness. In particular, under the condition L = (n + 1 2 )L D , where n is an integer number, the equilibrium β tb = 0 and gradually varies with the thickness: for L = (n + The modulation of the out-of-plane component of the magnetization of the surface spin spiral should lead to appearances of the stray fields above the sample, which in turn can be detected by Magnetic Force Microscopy (MFM) in films as well as in bulk crystals. However, with ordinary MFM technique it seems to be hard to distinguish the StSS from the ordinary helical state, which also produces a stray field around the sample. Moreover, MFM is not able to detect the relative orientation of the surface spin spirals on the opposite surfaces which is desired for experimental revealing of the result of our theoretical calculations.

Stacked spin spirals
The most promising experimental technique, which may allow to detect the StSS seems to be Lorentz TEM. Note, typically the Lorentz TEM images provide information only about the in-plain component of the magnetization. Figure 5 shows the in-plane component of the magnetization averaged over the film thickness. Two cases with stripe-and square-like pattern presented in figure 5(a) and (c) corresponds to the spin structure shown in figure 3(a) and (d), respectively.

Method
To construct the phase diagram shown in figure 2(a), we performed an energy minimization for each of the states and compared their energy density, see for instance figure 2(b). A finite difference approximation was used to convert the Hamiltonian (1) into a function, which arguments are the components of unit vectors defined on the nodes of finite size grid with unit cell size ∆x × ∆y × ∆z. Such calculations have been done for certain set of reduced thicknesses, L/L D , defined by the geometry of the simulated domain, material parameters and applied fields varying in the range 0 ≤ H/H D ≤ 1.5. For the function minimization we use a nonlinear conjugate gradient method implemented on NVIDIA CUDA architecture. In order to achieve the highest performance of the code and to keep the constraint n 2 = 1, we use the method of adaptive stereographic projections introduced in Ref. [19].
The size of the simulated domain along each of the spatial directions is defined as l x = N x ∆x, l y = N y ∆y, l z = (N z − 1)∆z, where N x , N y and N z are the number of nodes fixed to 128, 256 and 256...512, respectively (up to 16×10 6 nodes in total). For the fixed thickness, L, the value of ∆z is chosen such to satisfy the equality L = (N z − 1)∆z, while ∆x and ∆y are defined self-consistently to identify such l x and l y , which correspond to the lowest average energy density of the state. We assume periodic boundary conditions in the xy-plane and open boundary conditions along z-axis at z = 0 and z = L. Thereby, the procedure of energy minimization for each point on the phase diagram and for each of the states consist of an a priori unknown number of alternating steps: i) the direct energy minimization for given ∆x, ∆y, ∆z and ii) the small variation of ∆x and ∆y. In most general case, the variation of ∆x and ∆y is assumed to be independent. In particular, it is important for the proper energy minimization of StSS state, which is in general characterized by the modulations in both x-and y-directions, see figure 3. For the case of helical spin spiral and assuming k x-axis, the optimal ∆x minimizes the energy of the system when l x equals to the equilibrium period of the helix. Because of the homogeneity of the helical state in the direction perpendicular to k-vector, the variation in ∆y does not affect the energy density of the solution at all, and the problem reduces to a quasi-two-dimensional one.
The procedure based on direct minimization of the functional discussed above had been applied for relatively thin films, L ≤ 4L D . For the thick films, we take into account the exponential decay of surface induced modulation [21], which becomes negligibly small on the distance from the surface > L D . Thereby, the average energy density of any modulated state in the extended (infinite along the x-and y-directions) thick film with L > 4L D with a very high precision can be approximated as where e v is the average volume energy density for a particular state, which is assumed to be homogeneous along the applied field and e s is the average surface energy density of the twisted state. The In the range of thicknesses 4L D < L < 8L D we performed the calculations of the phase transitions with both approaches: i) direct energy minimization in whole volume of the layer and ii) according to Eq. (3). For this range of thicknesses we found identical results with both approaches.
In conclusion, one has to mentioned that the finding of the global energy minima for the functional (1) is by no means a straight-forward task due to the following reasons. Because of the complex energy landscape of 3D model with large energy barriers between the different states, the practically used algorithms of minimization provide only the local energy minimum close to the initial state, which in general may not correspond to the global minimum. To define the energetically dominant state one has to calculate and compare the energies of all competing phases. When the set of such tested states is incomplete it is impossible to identify the phase transitions properly. For instance when the earlier unknown StSS state is ignored the final phase diagram is incorrect [20].

Conclusions
We have presented the complete phase diagram for a film of isotropic chiral magnets, which allowed us to identify the range of existence of the equilibrium skyrmion lattice in applied magnetic fields and varying thicknesses. In particular, we found the critical thickness of the film above which the skyrmions never appear as the ground state of the system. We predict the existence of a new modulated state of stacked spin spirals, which occupies a wide range of the phase diagram for thin films as well as for bulk crystals. Such a state represents the coexistence of the conical spiral and surface induced spirals localized near the surface with finite penetration depth. We are confident that the presence of such states as well as the validity of the phase diagram itself can be confirmed with various experimental techniques.