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

Chiral magnets (ChMs) are a distinct class of magnetic crystals. Contrary to the classical ferro- and antiferromagnets the ground state of ChM is an incommensurate homochiral spin spiral—the 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_xGe [5], Mn${}_{1-x}$Fe_xSi [6], Fe${}_{1-x}$Co_xSi [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.

Zoom In Zoom Out Reset image size

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 [8, 9].

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 [10], 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 real bulk samples that are always covered by surfaces and thus are in the context of this paper semi-infinite solids. We provide a comprehensive description of this phase together with our suggestions how such a state can be experimentally detected.

It is worth mentioning that the existence of the stable magnetic skyrmions in magnetic crystals with DMI has been theoretically predicted by Bogdanov and Yablonskii in 1989 [11], while metastable solutions had been discussed right after in [12]. Nevertheless, it took more then twenty years to discover such particlelike states 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. [13], 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 [14]. However, the real nature of the A-phase still is under debate [14–16]. The theoretical models predict that in bulk crystals of ChMs, the magnetic skyrmions can be stabilized due to the strong cubic or uniaxial anisotropy [13, 17] or special crystal symmetry, which suppress the formation of the conical phase [11, 18]. 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 B20-alloys where the skyrmions are observed in a wide range of temperatures, much lower than the ordering temperature.

As has been shown in [19] the key to the understanding of the mechanism of skyrmion stabilization in thin films of ChM are presence of the free surfaces and the 3D 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 bobber—particlelike objects localized in all three spatial directions near the surface of the chiral magnet [20], see figure 1(f). Note, because of the presence of a singularity the chiral bobbers appear exclusively as a metastable state, at least at low temperatures. Contrary to topologically analogous states known in superfluid A-phase of He³ as vortices with the ends [24, 25], which are known to be unstable [26], the chiral bobbers in ChMs appear as a stable solution.

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.

The basic model for IChM include three main energy terms: the Heisenberg exchange interaction, the DMI and Zeeman energy term [27, 28]:

$$\begin{eqnarray}&&{ \mathcal E }={{ \mathcal E }}_{0}+{\displaystyle \int }_{V}{ \mathcal A }({{\partial }_{x}{\bf{n}}}^{2}+{{\partial }_{y}{\bf{n}}}^{2}+{{\partial }_{z}{\bf{n}}}^{2})+{ \mathcal D }\;{\bf{n}}\cdot [{\rm{\nabla }}\times {\bf{n}}]+{{HM}}_{{\rm{s}}}(1-{n}_{z}){\rm{d}}{\bf{r}},\end{eqnarray} \tag{ 1 }$$

where ${\bf{n}}\equiv {\bf{n}}({\bf{r}})$ is a continuous unit vector field defined everywhere except at the singular points. ${{ \mathcal E }}_{0}$ is the energy of the saturated ferromagnetic state. ${ \mathcal A }$ and ${ \mathcal 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 [29], and advanced models, which take into account the exact B20 crystal symmetry [30, 31]. Here we do not consider the contribution of dipole–dipole interaction, (i) because in the basic model of ChMs, this interaction plays the role of relatively small perturbation, while the main energy balance is defined by competition of the leading energy terms of Heisenberg exchange and DMI [8, 12], (ii) in order to keep consistency with earlier results obtained in frame of basic model [17–19].

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 [10, 18]:

$$\begin{eqnarray}&&{L}_{{\rm{D}}}=4\pi \displaystyle \frac{{ \mathcal A }}{{ \mathcal D }},\ {H}_{{\rm{D}}}=\displaystyle \frac{{{ \mathcal D }}^{2}}{2{M}_{{\rm{s}}}{ \mathcal A }}.\end{eqnarray} \tag{ 2 }$$

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.

3.1. 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 ${ \mathcal A },{ \mathcal D }$ and M_s, see equation (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.

Zoom In Zoom Out Reset image size

The solid lines in figure 2(a) correspond to the first-order phase transitions between helical spin spiral (red), skyrmion lattice (gray), 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\;\lt \;L/{L}_{{\rm{D}}}\;\leqslant \;4$, had been presented earlier in [20] . An attempt to generalize such a phase diagram for a wider range of thicknesses, $1\leqslant \;L/{L}_{{\rm{D}}}\lt \;\infty$, had been announced recently [21]. The phase diagram presented in [21] 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 triple point I defines the critical thickness of the film, ${L}^{*}=8.17{L}_{{\rm{D}}}$, above which the skyrmions may exist only as a metastable state. For instance, for MnSi (${L}_{{\rm{D}}}=18$ nm) and FeGe (${L}_{{\rm{D}}}=70$ nm) it gives ${L}_{\mathrm{MnSi}}^{*}\approx 150\;{\rm{nm}}$ and ${L}_{\mathrm{FeGe}}^{*}\approx 570\;{\rm{nm}}$. With decreasing thickness the range of existence for the skyrmion lattice in an applied magnetic field becomes wider. This fact reflects the relative contribution of the surface induced twists, which increases with decreasing thickness. There is another critical point, for $L/{L}_{{\rm{D}}}\lt 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 ${\bf{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 [22] and has been confirmed recently by directly observation with the Lorentz transmission electron microscopy (TEM) [23].

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}_{{\rm{D}}}\to 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}_{\mathrm{nv}},{D}_{2{\rm{d}}}$ and S₄ [11, 18].

3.2. Stacked spin spirals

A wide range of the phase diagram is occupied by the newly found StSS. The triple point II defines the limiting thickness above which the StSS may appear as the global energy minimum.

Figure 3 illustrates the complex spin structure of the StSS obtained by direct minimization of the functional (1). The StSS represents the coexistence of the conical phase in the bulk of the sample and the quasi-helical modulation of magnetization localized in the vicinity of the surface of the sample. Such free surface induced modulations have finite penetration depth and appear on both the top and bottom surfaces. It exhibits a mixed helical- and cycloidal-like modulation as shown in figure 4(a). The clock-wise rotation of magnetization in the helical-like part is chosen as direction of the wave vector ${\bf{k}}$ for such a complex spin spiral, see figure 4(c).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We found that the relative orientation of the propagation vectors of the spin spirals on the top and bottom surfaces, ${{\bf{k}}}_{{\rm{t}}}$ and ${{\bf{k}}}_{{\rm{b}}}$, respectively, is thickness dependent. In other words, the angle ${\beta }_{\mathrm{tb}}$ between ${{\bf{k}}}_{{\rm{t}}}$ and ${{\bf{k}}}_{{\rm{b}}}$ for the equilibrium StSS varies as function of the film thickness. In particular, under the condition $L=\left(n+\tfrac{1}{2}\right){L}_{{\rm{D}}}$, where n is an integer number, the equilibrium ${\beta }_{\mathrm{tb}}=0$ and gradually varies with the thickness: for $L=\left(n+\tfrac{3}{4}\right){L}_{{\rm{D}}},{\beta }_{\mathrm{tb}}=90^\circ$, for $L=(n+1){L}_{{\rm{D}}},{\beta }_{\mathrm{tb}}=180^\circ$, and for $L=\left(n+1\tfrac{1}{4}\right){L}_{{\rm{D}}},{\beta }_{\mathrm{tb}}=90^\circ$. Figure 3(a) illustrates the case of ${\beta }_{\mathrm{tb}}=180^\circ ,L=4{L}_{{\rm{D}}}$, while (d) illustrates the case of ${\beta }_{\mathrm{tb}}=90^\circ ,L=4\tfrac{1}{4}{L}_{{\rm{D}}}$. Note, because of the relatively weak energy dependence upon the ${\beta }_{\mathrm{tb}}$ in the case of thick films, one should expect a significant influence of the magnetocrystalline anisotropy to the orientation of the vectors ${{\bf{k}}}_{{\rm{t}}}$ and ${{\bf{k}}}_{{\rm{b}}}$ with respect to certain crystallographic directions. In case of very thin films, $L/{L}_{{\rm{D}}}\lesssim 3$ where StSS state becomes metastable and surface modulations significantly interfere each other the equilibrium orientation of the spirals in opposite surfaces does not satisfy above rules.

Figure 5(a) illustrates the in-field dependence of the period of the surface induced modulations, yellow line, in comparison to the period of the helical and conical spirals, red and blue lines, respectively. These dependencies corresponds to the limiting case of the semi-infinite crystal $L/{L}_{{\rm{D}}}\to \infty$, when the interaction between opposite surfaces can be neglected, and modulations at only one surface have to be taken into account. Note, the period of surface induced spirals is always larger than that of ordinary helical spirals, ${P}_{\mathrm{StSS}}\geqslant 1.41{P}_{\mathrm{Hel}}$, and increases with the applied field. The dependence of the energy density for StSS as a function of the distance from the surface for different applied fields is shown in figure 5(b). It illustrates the exponential decay of surface induced modulations [22], which for $z^{\prime} \gt 1.5{L}_{{\rm{D}}}$, deep inside of the crystal, converges to the pure conical phase. Because of such a strong localization of the surface induced modulations for relatively thick films with $L\gtrsim 4{L}_{{\rm{D}}}$, the energy of the surface modulations becomes almost thickness independentthe modulations on opposite surfaces do almost not affect each other. For $H\lt 0.42{H}_{{\rm{D}}}$ and $L\gtrsim 4{L}_{{\rm{D}}}$ the energy of the StSS is always lower than the energy of the conical phase, see e.g. yellow line in figure 2(b). Therefore, such surface modulations should appear even near the surface of the bulk crystals when the applied field is pointing along the surface normal.

Zoom In Zoom Out Reset image size

We are fully aware that the basic model considered in the present work may not be met by the surfaces of real samples, which may exhibit effects like surface alloying, structural modifications, the formation of terraces and others. All this will change the electronic structure at the vicinity of the surface and thus the parameters entering our model. Even for an ideal surface the parameters entering the model will be changed within a few atomic layers in the close vicinity to the free boundary since the electronic structure of such surface atomic layers differ from the electronic structure of internal atoms. Because the coupling constants for the exchange interaction and spin–orbit interaction are changed, we may expect in general a small modification of the absolute values of phase transition lines and critical points on the phase diagram. However, the contribution of this kind of surface effect in first approximation should be proportional to the ratio between thickness of such surface layer λ on the electronic scale and the penetration depth of surface induced modulations, $\lambda /{L}_{{\rm{D}}}$. Typically λ is in order of maximal 10 atomic lattice constants a, while L_D is in order of few tens or hundreds of a. As a result the contribution of the surface atomic layers should be negligibly small with respect to the chiral surface twist effect for the most of ChMs with long periodical modulations and our finding is a generic feature of such systems.

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. Note, the absolute values of periodicity for the helical spiral and surface modulations, nevertheless, allow distinguishing between these two phases, see figure 5(a). 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.

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-plane component of the magnetization. Figure 6 shows the in-plane component of the magnetization averaged over the film thickness. Two cases with stripe- and square-like pattern presented in figures 6(a) and (c) corresponds to the spin structure shown in figures 3(a) and (d), respectively.

Zoom In Zoom Out Reset image size

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 ${\rm{\Delta }}x\times {\rm{\Delta }}y\times {\rm{\Delta }}z$. Such calculations have been done for certain set of reduced thicknesses, $L/{L}_{{\rm{D}}}$, defined by the geometry of the simulated domain, material parameters and applied fields varying in the range $0\leqslant H/{H}_{{\rm{D}}}\leqslant 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 ${{\bf{n}}}^{2}=1$, we use the method of adaptive stereographic projections introduced in [20].

The size of the simulated domain along each of the spatial directions is defined as ${l}_{x}={N}_{x}{\rm{\Delta }}x,$ ${l}_{y}={N}_{y}{\rm{\Delta }}y,{l}_{z}=({N}_{z}-1){\rm{\Delta }}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\times {10}^{6}$ nodes in total). For the fixed thickness, L, the value of ${\rm{\Delta }}z$ is chosen such to satisfy the equality $L=({N}_{z}-1){\rm{\Delta }}z$, while ${\rm{\Delta }}x$ and ${\rm{\Delta }}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 ${\rm{\Delta }}x,{\rm{\Delta }}y,{\rm{\Delta }}z$ and (ii) the small variation of ${\rm{\Delta }}x$ and ${\rm{\Delta }}y$. In most general case, the variation of ${\rm{\Delta }}x$ and ${\rm{\Delta }}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 ${\bf{k}}\parallel x$-axis, the optimal ${\rm{\Delta }}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 ${\bf{k}}$-vector, the variation in ${\rm{\Delta }}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\leqslant 4{L}_{{\rm{D}}}$. For the thick films, we take into account the exponential decay of surface induced modulation [22], which becomes negligibly small on the distance from the surface $\gtrsim {L}_{{\rm{D}}}$. Thereby, the average energy density of any modulated state in the extended (infinite along the x- and y-directions) thick film with $L\gt 4{L}_{{\rm{D}}}$ with a very high precision can be approximated as

$$\begin{eqnarray}&&e(H,L)={e}_{{\rm{v}}}(H)+\displaystyle \frac{2{e}_{{\rm{s}}}(H)}{L},\end{eqnarray} \tag{ 3 }$$

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. For instance for StSS state e_s equals to the total area shown as red and blue domains in figure 5(b). Note, red domains correspond to the energy gain and have negative sign while blue domains correspond to the energy loss and have positive sign. The ${e}_{{\rm{v}}}={{ \mathcal E }}_{{\rm{v}}}/V$, where ${{ \mathcal E }}_{{\rm{v}}}$ is the total energy calculated on the domain with $L=4{L}_{{\rm{D}}}$ in the absence of free surfaces meaning periodic boundary conditions in all three spatial directions; V is the volume of the simulated domain. Then, the total energy ${{ \mathcal E }}_{{\rm{v}}}^{\prime }$ has been calculated on the same domain with open boundary condition at z = 0 and z = L. The surface energy density e_s has been identified from the difference ${{ \mathcal E }}_{{\rm{v}}}-{{ \mathcal E }}_{{\rm{v}}}^{\prime }$. For each phase with surface induced modulations, the value of e_v and e_s has to be calculated ones for each fixed value of H and then the energy density in whole range of thicknesses $L\geqslant 4{L}_{{\rm{D}}}$ can be found according to equation (3).

In the range of thicknesses $4{L}_{{\rm{D}}}\lt L\lt 8{L}_{{\rm{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 equation (3). For this range of thicknesses we found identical results with both approaches.

In conclusion, one has to mention 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 [21].

We have presented the complete phase diagram for a film of isotropic ChMs, 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.

The research of FNR and ABB was carried out within the state assignment of FASO of Russia (theme Quantum No. 01201463332). The work of FNR was supported in part by RFBR (project No. 14-02-31012).