Search tips
Search criteria 


Logo of scirepAboutEditorial BoardFor AuthorsScientific Reports
Sci Rep. 2016; 6: 29126.
Published online 2016 July 5. doi:  10.1038/srep29126
PMCID: PMC4932608

Exotic skyrmion crystals in chiral magnets with compass anisotropy


The compass-type anisotropy appears naturally in diverse physical contexts with strong spin-orbit coupling (SOC) such as transition metal oxides and cold atomic gases etc, and it has been receiving substantial attention. Motivated by recent studies and particularly recent experimental observations on helimagnet MnGe, we investigate the critical roles of this compass-type anisotropy in modulating various spin textures of chiral magnets with strong SOC, by Monte Carlo simulations based on a classical Heisenberg spin model with Dzyaloshinsky-Moriya interaction and compass anisotropy. A phase diagram with emergent spin orders in the space of compass anisotropy and out-of-plane magnetic field is presented. In this phase diagram, we propose that a hybrid super-crystal structure consisting of alternating half-skyrmion and half-anti-skyrmion is the possible zero-field ground state of MnGe. The simulated evolution of the spin structure driven by magnetic field is in good accordance with experimental observations on MnGe. Therefore, this Heisenberg spin model successfully captures the main physics responsible for the magnetic structures in MnGe, and the present work may also be instructive to research on the magnetic states in other systems with strong SOC.

Spin-orbit coupling (SOC) plays an important role in condensed matter physics, especially in recently addressed quantum spin Hall effects and related topological orders1,2. So far, magnetic materials with lacking inversion symmetry in B20 cubic crystal structure provide an alternative arena for studying the emergent phenomena in magnets of strong SOC. Prominent members encompass metallic compounds like MnSi3,4,5 and FeGe6, semiconductor Fe1–xCoxSi7,8, and multiferroic insulator Cu2OSeO39. In these materials, the competitions between strong ferromagnetic (FM) exchanges, weak Dzyaloshinsky-Moriya (DM) interactions, and additional magnetic anisotropies originating from high-order interactions, generate a long-wavelength helical spin order which propagates along some specific directions favored by these magnetic anisotropies3,4,5.

The roles of magnetic anisotropies are usually delicate but substantial in determining the magnetic ground states in some cases. Generally, these anisotropy terms are the weakest in energy scale among the above mentioned interactions, and the favored directions of helical spin orders are materials-dependent. For example, the helical order in MnSi observed in neutron scattering experiments is characterized by a set of Bragg peaks situated on a sphere surface, suggesting that the spin alignment of the helical state is locked along the cubic space diagonal <111> directions, due to the effect of high-order crystal anisotropies3,4,5. On the other hand, a perpendicular external magnetic field may transform the single twisted helice to the energetically favorable hexagonal skyrmion crystal (SkX) structure4,8. This intriguing topological structure has been receiving tremendous interest due to its unusual magnetic and transport properties, especially the topological Hall effect10,11,12.

Nevertheless, recent experiments on B20-type cubic MnGe11,13 revealed the field dependence of the topological Hall effect which is distinct from those observed in other B20-type compounds, suggesting a zero-field ground state skyrmion lattice very different from often observed hexagonal SkX. So far, no recognized theory for this unusual skyrmion lattice is available. In addition, the helical wavelength λ changes from 3 nm to 6 nm with temperature, which is much smaller than that in conventional B20-type compounds with typical λ of 18–90 nm. Therefore, a different mechanism for the spin structure in this generic B20 MnGe is expected. One possible interpretation of the variation in small-angle neutron scattering (SANS) patterns with magnetic field is the presence of an easy axis of magnetization along the <100> direction13. This implies that MnGe may have strong SOC and consequent non-negligible and high-order magnetic anisotropy due to SOC, so that the modulation vector q is locked along the <100> direction. In this sense, the appearance of exotic skyrmion state is easy to understand.

Another scenario for the existence of high-order anisotropy terms in B20-type transition-metal silicides and germanides is associated with the weak itinerant-electron ferromagnetism14. The outermost electrons in these compounds can be described by an extended Hubbard model on the two-dimensional (2D) lattice. The Hamiltonian consists of the nearest-neighbor natural hopping (t0), SOC-induced hopping (tSO), and Hubbard repulsion (U)15,16,17,18:

An external file that holds a picture, illustration, etc.
Object name is srep29126-m1.jpg

where An external file that holds a picture, illustration, etc.
Object name is srep29126-m2.jpg is the operator creating a spin-σ (σ = ↑, ↓) electron at site i, and An external file that holds a picture, illustration, etc.
Object name is srep29126-m3.jpg denotes the elements of the Pauli matrix with superscript η = x, y. In the large U condition, one can obtain an effective spin Hamiltonian, by a combination of the Heisenberg exchange (J), DM interaction (D), and quantum compass anisotropy (Ac)15,16,17,18,19:

An external file that holds a picture, illustration, etc.
Object name is srep29126-m4.jpg

where An external file that holds a picture, illustration, etc.
Object name is srep29126-m5.jpg is the spin vector with μ [equivalent] (x, y, z). By defining t = (t02+tSO2)1/2, tanθ = tSO/t0, and J0 = 4t2/U, one can obtain J = J0 · cos2θ, D = J0 · sin2θ, and Ac = J0 · (1  cos2θ)16. For magnetic systems with strong SOC, an easy-plane anisotropy term, i.e. the so-called “compass anisotropy (Ac > 0)”, appears naturally15,16. This Hamiltonian has recently been used to deal with ultra-cold atoms with artificial SOC on a 2D optical lattice or three-dimensional (3D) counterpart17,18.

In fact, in the context of solid state physics, hybrid models interpolating between Heisenberg models and compass models are often proposed20,21,22,23. Such hybrid models resulting from Hubbard models are relevant to describe the superexchange interactions in transition metal systems with strong SOC, such as those containing 4d and 5d ions Rh, Ru, Os, and Ir, where the effective moments carry both the orbital and spin characters20,22,23. An example can be found in recent experiments demonstrating the coexistence of superconductivity and ferromagnetism on the interfaces between LaAlO3 and SrTiO3 (2D electron gas)24. Interestingly, relevant magneto-transport studies show a strong and gate-tunable Rashba SOC for these 2D conduction electrons, arising from the broken inversion symmetry at the interface25. Subsequently, the microscopic mechanism for the interfacial magnetism using a microscopic model (derived from an extended Hubbard model with Rashba SOC) was investigated15,16. Both the numerical and analytical calculations predict a long-wavelength spiral ground state with a SOC-dependent pitch under zero magnetic field, while hedgehog-like skyrmions appear over a broad range of magnetic field. While the relevant mechanisms for these emergent phenomena have been disputed, the possible role of the compass anisotropy thus becomes interested to some extent15,16,17,18,20,22.

From the viewpoint of energy scale, the compass anisotropy is usually ignored in the literature, since it is an even higher order than the DM term which is already insignificant. We discuss this exceptional issue from several levels:

(1) First, considering a chiral magnet with commonly weak SOC (tSO/t0 [double less-than sign] 1), the contribution of compass anisotropy term to the total energy is on an order of magnitude of Ac ~ J · (tSO/t0)2 which is much smaller than the DM term D ~ J · (tSO/t0) [double less-than sign] J16. However, MnGe has strong SOC and this compass anisotropy term can no longer be trivial and must be taken into account, for its contribution to the energy is comparable to the DM term (~D2/J)15,16.

(2) Second, the effect of uniaxial anisotropy As(Siz)2 on the helical state and SkX in B20-type helimagnets was once studied26. In phenomenological sense, this uniaxial anisotropy arises from the single-ion or dipolar shape anisotropy, and thus can be either an easy-axis (As < 0) or hard-axis (As > 0) anisotropy27. To this stage, the effective anisotropy in these materials is governed by A = Ac + As. In some cases, the uniaxial anisotropy is vanished or much smaller than the compass term, which is typically true in MnGe and those magnetic systems with strong SOC mentioned above15,16. This implies that the compass anisotropy can be a dominant ingredient in the effective anisotropy term in these cases, which has not been well recognized earlier11,13,15,16.

(3) Third, it is noted that the anisotropy ingredients would possibly be sensitive to external stimuli. While the uniaxial anisotropy can be tuned by strain27, the compass anisotropy may be modulated largely too by external stimuli (e.g. external electric field), as motivated by the tunable SOC in 2D electron gas systems28,29. These allow possibilities to explore the spin states in chiral magnets with tunable A (either uniaxial or compass or both). Obviously, along the line supported by the above three-level discussion, one has sufficient reasons for a theoretical study on the helimagnets with compass anisotropy like MnGe.

In this work, we investigate the spin ground states in such chiral magnets with multifold interactions, using Monte Claro (MC) simulation30. Specific attention is paid to the roles of the compass anisotropy. We first explore the phase diagram in the (Hz, A)-space, which demonstrates the roles of the compass anisotropy in modulating the spin states. Then, we suggest an intriguing alternating half-skyrmion (HSk) and half-anti-skyrmion (HASk) crystal structure (i.e., a state in the phase diagram) to be a candidate for the zero-field ground state of MnGe compound, confirmed by the quite good qualitative consistence of the MC simulations with experiments. We believe that this study provides a theoretical guide to understand the magnetic structures and their evolutions in those helimagnets with strong SOC.


Model and Simulation

Considering the helical and skyrmion spin structures in helimagnets MnGe, MnSi, Fe1−xCoxSi etc, it does not lose a generality to start from a 2D L × L square lattice with periodic boundary conditions. One Heisenberg spin Si is imposed on each site i. The Heisenberg exchange interactions, DM interactions, compass anisotropy terms, and Zeeman term are considered. The Hamiltonian Eq. (2) is thus re-written as:

An external file that holds a picture, illustration, etc.
Object name is srep29126-m6.jpg

where field H is along the z-axis or –x-axis or stated elsewhere. The uniaxial anisotropy favors the hard-axis along the z-axis at As > 0, or the easy-axis along the z-axis at As < 0. Here, Asis a minor quantity in comparison with Ac. In this sense, one may reckon that the synthetic anisotropy term (A) still prefers the in-planar 90° compass-type symmetry. We will simply treat this synthetic easy-plane anisotropy as the compass anisotropy hereafter.

In the classical approximation, we treat spins Si as classical vectors and aim at minimizing the energy for finding the ground-state spin configurations {Si}. It is noted that the classical MC simulation has been used to explore the phase diagrams of effective spin models in the similar context, e.g. chiral magnets8,31, and even ultra-cold atomic systems17,18. Therefore, we restrict our discussion on the classical approximations, while a full quantum treatment would be needed to have a deeper understanding of quantum spin models, which is beyond the present work.

For simplicity, we use dimensionless units for convenience: spatial length is scaled by the lattice spacing a. Parameters D, A, and H are quantified in the units of J/a, J, J · S, respectively. The spin moment S [equivalent] 1 and a = 1 are taken without loss of generality. In this case, the spin phase texture is determined by ratios between D, A, and H. In the simulation, we choose ratio D/J = √6, yielding the helical wavelength λ ~ 8a for the spin structure in MnGe. Note that the typical helical wavelength λ is ~4.0 nm for MnGe, and the lattice constant a is ~0.48 nm32, which is on the same order of magnitude as that of the helical states obtained in our simulations. The simulation algorithm and details of the calculation procedure are presented in the Method section below.

The as-simulated spin structures are characterized by several physical quantities, as suggested in literature31. First, the spatial symmetry of spin structure is reflected in the Bragg intensity pattern |S(q)|2, with reciprocal momentum S(q) obtained from the Fourier transformation of the spin configuration:

An external file that holds a picture, illustration, etc.
Object name is srep29126-m7.jpg

where N = L2 is the total number of spins and ri is the spatial coodinate.

To characterize swirling structure of the skyrmions lying on the xy-plane, we introduce the local density of skyrmions χi at lattice site i, defined as:

An external file that holds a picture, illustration, etc.
Object name is srep29126-m8.jpg

which is the discretization counterpart of the skyrmion density S·([partial differential]xS × [partial differential]yS)/4π for the continuum case4. The summation of χi over the extended coordinate space gives the so-called topological winding number (skyrmion number) as χ = iχi, which is proportional to the topological Hall resistivity An external file that holds a picture, illustration, etc.
Object name is srep29126-m9.jpg . One single skyrmion cell contributes one unit to χ in the continuum limit.

The helicity of a skyrmion state or a helical state is closely related to the spin chirality via the relativistic SOC. We define the helicity as33:

An external file that holds a picture, illustration, etc.
Object name is srep29126-m10.jpg

where the sign and magnitude of γ reflect respectively the spin swirling direction (i.e., left or right handedness) and degree of swirl in a proper screw spin structure.

Moreover, we define the out-of-plane magnetization, in-plane magnetization, and total magnetization, respectively as:

An external file that holds a picture, illustration, etc.
Object name is srep29126-m11.jpg

where <…> refers to the thermal and configuration averaging at a given temperature T.

Low temperature phase diagram

Extensive simulations over a broad region on the (A, Hz) plane generate a low temperature phase diagram, as summarized in Fig. 1. Here Hz denotes the field H along the z-axis. Some typical spin structures of these phases are presented in Fig. 2. To illustrate the spin configurations, we use the color map to scale the spin components along the z-axis (out-of-plane) Siz, and use the spin vectors (arrows) to describe the on-plane xy components Sixy.

Figure 1
Low temperature (T = 0.01) phase diagram of the spin model in Eq. (3) with D = √6J, and H along the z-axis (Hz).
Figure 2
Plot of some typical spin configurations arising in the phase diagram (Fig. 1): (a) HP state at (A, Hz) = (0.2, 0.0), (b) SkX phase at (A, Hz) = (0.2, 2.6), (c) SC1 phase at (A, Hz) = (4.2, 0.0) and (d) ...

For very low Hz ~ 0 in Fig. 1, the result shows the spin structure transition from helical phases (HP) to spin crystal (SC) phases, as the anisotropy A increases from zero. Upon applying of an intermediate Hz, the HP and SC1 phase evolve respectively into the SkX phase and SC2 phase, as shown in Fig. 2(b,d). At high Hz, the lattice gradually gives away to the FM phase, with all spins orientating along the z-axis. This phase transition sequence is consistent with the one in previous studies31.

In Fig. 2(c,d), one may find from the Bragg patterns that these crystal structures are constructed primarily from the superimposition of two helices, one oriented along the x-axis and the other along the y-axis, as characterized by two pairs of Bragg peaks. Here, the equality of these two pairs of Bragg intensities can be used as a criterion to classify the chiral crystal states: the SC1 phase (unequal intensities) and SC2 phase (equal intensities), following earlier studies31. Therefore, the increasing compass anisotropy leads to the transformation from single-q spin structures (i.e. HP) to double-q spin structures (i.e. SC phases). In the following of this work, we emphasize in the SC phases and investigate their topological properties and magnetic-field-induced evolution behaviors.

The SC2 phase: candidate for the ground state of MnGe

Now we investigate the spin crystal phases, i.e. the SC1 and SC2 phases shown in Fig. 2(c,d). These spin crystal phases were proposed theoretically34, and they have been confirmed by the recent experimental observations on MnGe11,13. The neutron scattering revealed the multiple-q helix structure in MnGe in response to increasing H, with the spin helix q vectors locked along the <100> direction due to an additional anisotropy originating from the SOC. More importantly, given a square or cubic SkX structure, the Berry phase calculations fit quite well the observed H-dependent topological Hall resistivity. In contrast, if a hexagonal SkX structure in MnSi and Fe1–xCoxSi is assumed, the calculations fail to fit the observed topological Hall resistivity as a function of H. Therefore, the MnGe most likely prefers the square or cubic SkX structures rather than a hexagonal SkX structure13.

The square SkX structure or the 2D counterpart of cubic SkX structure is just the SC2 phase presented in Fig. 2(d). However, knowledge on microscopic mechanism for the possible SC2 phase in MnGe is still lacking13. The present work proposes that it is just the compass anisotropy originating from the SOC responsible for formation of the SC phases, noting again that MnGe has strong SOC. To proceed, we start from a SC structure on a 2D plane, which can be viewed as a superimposition of two single helical states34:

An external file that holds a picture, illustration, etc.
Object name is srep29126-m12.jpg

where I [equivalent] I(x, y) is a normalized factor, B (0 < B  1) is a variable associated with the compass anisotropy term A, which stands for the compression of the helical ordering along the z-axis brought about by the easy-plane compass anisotropy term. C1 and C2 represent the specific weight of the two helices, as characterized by the relative intensities of the two pairs of the Bragg spots shown in Fig. 2(c,d).

The SC1 and SC2 phases as the ground states at H~0

We first study the spin structures and possible topological properties of the SC1 and SC2 phases under H = 0, corresponding to the Eq. (8) with C1  C2 and C1 = C2, respectively. The typical spin configurations of the two phases are plotted in Fig. 3(a,b) respectively. While the SC1 phase is relatively trivial, we concentrate on the SC2 phase. We find a periodic array of “nodes” (magnetization nodes) in the SC2 phase, where the spin vector Si = 035. These singular points can be easily seen from Eq. (8) if one sets (sinqx = 1, cosqy = −1) or (sinqx = −1, cosqy = 1). The hard-spin constraint |Si| = 1 are used in the simulation, where the singularities are naturally forbidden. In fact, however, it is possible that the spins in real chiral magnets are soft due to their interactions with the metallic host or the averaging over fast fluctuations, leading to the possible points of vanishing magnetization. The electronic structure of this nodal topological lattice had been studied recently36. However, we find that the magnetic nodes don’t impose substantial impact on the field dependent behaviors and topological properties of the SC2 phase, and here we no longer give more discussion on the nodes. The SC2 phase consists of alternative alignment of two square flux units, with two specific units enclosed by the dash lines in Fig. 3(b). We calculate the topological winding number χ of such two units by adopting a numerical procedure. First, the 2D lattice is discretized on a square mesh with grids up to N × N, and then the χi at each grid point is calculated using Eq. (5). Second, all the χi over the defined region are summed, where an enough big N is taken so that the error is on the order of An external file that holds a picture, illustration, etc.
Object name is srep29126-m13.jpg 37.

Figure 3
Plots of (a) SC1 structure given by the expression Eq. (8) with specific parameters q = 2π/30, C1 = 0.8, C2 = 0.2, and B = 0.9; and (b) SC2 structure with B = 0.5, ...

The numerical results show χ = 1/2 for one flux unit, and χ = −1/2 for another one. Note that the calculated results don’t depend on factor B within 0.0 < B  1.0. The topological winding number χ = 1/2 or −1/2 for the square flux indicate that this unit is a half-anti-skyrmion (HASk) structure or a half-skyrmion (HSk) structure, with the schematic spin configurations shown in Fig. 3(c,d), respectively. Therefore, the SC2 phase at H ~ 0 is a spin crystal structure consisting of alternating HSk and HASk units.

It is known that a skyrmion can be topologically mapped into a hairy sphere (or a hedgehog). The vortex-like skyrmion texture found in MnSi and Fe1–xCoxSi is topologically equivalent to such a hairy sphere38. Applying this topology property to the SC2 phase here, one may unfold the HSk and HASk units and map them respectively onto the upper hemisphere and lower hemisphere of the hairy sphere, as shown in Fig. 3(c)~(f). In the other words, the HSk and HASk units can be obtained from splitting a skyrmion into two halves. In this sense, the compass anisotropy causes the fractionalization of the skyrmions.

Evolution of the SC2 phase in response to Hz

Subsequently, we investigate the evolution of the SC2 phase in response to the cycling of H. To compare with experiments11,13, we focus on the cases of H along the z-aixs (Hz) and –x-axis (H−x). All the simulations start from the initial lattice which is equilibrated at H = 0, D = √6, A = 6.0, and T = 0.01, corresponding to the zero-field cooling in experiments. A series of parameters characterizing the SC2 phase in response to increasing Hz are plotted in Fig. 4. Here the Hz is varying gradually from zero to a value big enough for saturating the FM phase and then back to zero.

Figure 4
Plots of various physical quantities as a function of Hz: (a) <χ>, (b) <γ>, (c) thermal-averaged energy per spin <E>, (d) <Mxy>, (e) <Mz>, and (f) < ...

We first look at the evolutions in the Hz-increasing half loop. Figure 4(a) shows that the thermal-averaged <χ> gradually drops from zero to a negative maximal of −5.0 at Hz = 5.5, implying a broken balance of the χ-contributions from the HSk and HASk units in the SC2 phase. In the simulation, it is found that the positive skyrmion density in the HSk units reduces markedly, leading to a negative net spin chirality, with the spin structure at Hz = 5.5 shown in Fig. 4(g). Further increasing of Hz drives the <χ> back and eventually to die away at Hz ~ 12.0 (corresponding to the FM state). The evolution from the SC2 phase to the FM state is not always continuous, but featured with a sharp jump in the <χ>-Hz curve at Hz ~ 10.3, accompanied with an abrupt expansion of the HASk unit size and also abrupt jump of the SC periodicity. The spin structures right below and above this jump point are presented in Fig. 4(h,i). Similar phenomenon was also observed in the melting of hexagonal SkX structure. In this process, the increasing Hz destructs the thermal-averaged helicity <γ> of the SC2 phase, as seen clearly by the <γ>-Hz curves shown in Fig. 4(b).

Then we look at the Hz-decreasing half loop. Different from the Hz-increasing half loop, the lattice returns back to a spin structure aligned along the <10> direction rather than the SC2 phase, as shown in Fig. 4(j). This spin structure is characterized by a zigzag pattern with striped domain of alternating Mz. The simulations on an extended 64 × 64 lattice confirm the existence of the zigzag striped domains. It is noted that no thermal hysteresis effect appearing in the evolution of magnetic state during the process of increasing and decreasing the magnetic field, as shown by the simulated <E>-Hz loop in Fig. 4(c). Therefore, the different magnetic states emerging during this evolution sequence suggest that these states are basically degenerate. Typically, the zigzag striped domains and the SC2 phase at Hz = 0 have the same energy. The zigzag spin pattern may align along one of <10> directions on the xy-plane at random. Considering the zigzag-patterned spin structure which aligns along the y-axis shown in Fig. 4(j), the average x component <Mx> of magnetization is ~0, leaving the nonzero <My>. However, the values of <Mx> and <My> reverse for the zigzag pattern aligned along the x-axis. Therefore, we use in-plane magnetization <Mxy> rather than <Mx> and <Mz> to characterize the formation of zigzag-patterned spin structure upon decreasing Hz. As shown in Fig. 4(d)~(f), the simulated <Mxy>-Hz, <Mz>-Hz, and <Mxyz>-Hz curves all confirm that the zigzag spin structure with magnetization on the xy-plane develops gradually upon decreasing Hz, resulting in the nonzero total magnetization <Mxyz> of this spin pattern.

The above simulated results find qualitative consistence with recent experimental observations. First, the simulated <χ>-Hz loop is in perfect agreement with measured An external file that holds a picture, illustration, etc.
Object name is srep29126-m14.jpg -Hz loop on polycrystalline MnGe at low temperature (5 K), noting that the topological Hall resistivity An external file that holds a picture, illustration, etc.
Object name is srep29126-m15.jpg is proportional to χ (see Fig. 2(d) in [11]). In addition, the measured <M>-Hz loop does show a linear behavior (see Fig. 4(a) in [11]), consistent with the simulated <Mz>-Hz dependence. Considering the random distribution of zigzag-patterned spin structure or SC2 domains in the xy-plane for polycrystalline MnGe, one may compare the simulated <Mz>-Hz curves with the measured <Mxyz>-Hz curves. Indeed, they coincide with each other quite well.

The Hz-driven evolution from the SC2 phase to the FM phase can be qualitatively explained, according to a theoretical formula34:

An external file that holds a picture, illustration, etc.
Object name is srep29126-m16.jpg

where M1 denotes the net magnetization induced by Hz and M2 represents the magnetization of the SC2 phase, factor B (0 < B  1) again stands for the compression of the helical ordering along the z-axis induced by the easy-plane compass anisotropy.

To proceed, we separate one HSk-HASk pair into two independent primitive cells. Then, we discretize each of the HSk unit and HASk unit on a 2D mesh with up to 200 × 200 grids for calculating their topological winding number χHSk and χHASk, as the χ-Mz dependences plotted in Fig. 5(a). Although the field Hz is not contained explicitly in Eq. (9), we may use the χ-Mz dependence to characterize the χ-Hz dependence, by considering that Hz can be scaled by Mz as Mz increases along with Hz. We can see that the χHSk drops steeply from 0.5 to 0.0 and the χHASk from −0.5 to −1.0 with increasing Mz from 0.0 to ~0.2, implying a sharp transformation from the HASk to an anti-skyrmion with χ = −1 and from the HSk to a flux (vortex) state with χ = 0. Consequently, a nonzero net spin chirality over the whole xy-plane emerges, as the typical spin structure at Mz = 0.55 shown in Fig. 5(b). Further increasing Mz from ~0.2 drives the SC2 phase to deform continuously and eventually the anti-skyrmion state is melt, resulting in χHASk ~ 0.0 at Mz ~ 0.65. The spin chirality becomes zero, corresponding to the FM phase.

Figure 5
(a) A plot of the χ-Mz dependence for characterizing the evolution of the SC2 phase in response to Hz, according to Eq. (9) with B = 0.5 chosen here. χ is counted in the region of a HSk (HASk) structure, denoted as χ ...

To this stage, our calculations suggest the evolution sequence from the SC2 phase to the FM phase upon increasing Hz as: (1) HASk unit  anti-skyrmion  flux state/vortex with zero net spin chirality  FM state; (2) HSk unit  flux state/vortex with zero net spin chirality  FM state. Likewise, when H is applied along the –z-axis, a transition from a HSk unit to a skyrmion is expected.

Evolution of the SC2 phase in response to H−x

As a complimentary to the last section, we also try to uncover the evolution of the SC2 phase in response to an in-plane magnetic field, e.g. along the –x-axis (H−x). The main simulated results are summarized in Fig. 6. Here, we perform the MC simulation on a 64 × 64 lattice. Such a big size of lattice makes it better to distinguish the spin structures, and especially the periodicity of these structures. In Fig. 6(a), the simulated <Mxy>-Hx dependence in the increasing Hx sequence exhibits the multi-step magnetization. The alternating HSk and HASk crystal structure gradually decays into the FM phase with spin aligned along the –x-axis, as some typical spin lattices shown in Fig. 6(d)~(h). The periodicity of the spin patterns can be measured from the Bragg intensity patterns in the upper insets, which are associated with the experimental SANS patterns13.

Figure 6
Plots of various physical quantities as a function of Hx: (a) <Mxy>, (b) <γ>, and (c) energy per spin <E>. In this process, spin configuration evolves from the initial SC2 state at H = 0 ...

In fact, every magnetization step corresponds to a specific spin pattern. At the low-Hx case, shown in Fig. 6(d), the SC2 structure is slightly deformed from the C4 symmetry, reflected in the two pairs of Bragg spots with different brightness. It should be noted that the experimental ring-shaped SANS pattern, formed in the polycrystalline powder, is generated by the randomly oriented multiple-q structures with the same magnitude of magnetic modulation vector (see Fig. 3(a) in [13]). Therefore, the ring pattern is associated with the formation of a periodically modulated magnetic structure (i.e., the SC2 phase here). Further increasing Hx triggers some stripe-like domains along the x-axis shown in Fig. 6(e) at Hx = 0.35. These stripes are mixture of the alternating HSk and HASk crystal structure and the FM state aligned along the –x-axis with inharmonic periodicity. The corresponding Bragg intensity pattern shows a pair of Bragg spots in the qx-axis and two pairs of Bragg spots in the qy-axis, in good agreement with the crescent-shaped SANS pattern observed experimentally (Fig. 3(b) in [13]). At very high Hx, the Bragg spots shrink and the lattice gradually gives away to the FM phase, also consistent with experiments13. Similar features can be found in the <γ>-Hx curve (Fig. 6(b)). When Hx is gradually ramped down back to zero, a spiral state emerges at Hx ~ 0.3. As shown in Fig. 6(i), the spiral propagates along the –x-axis with a good periodicity, as the variation in the Bragg patterns again shows good consistence with the variation of SANS patterns (see Fig. 3(e)~(h) in [13]).

It is noted from the <E>-Hx curves in Fig. 6(c) that, the spiral state is degenerate with the alternating HSk and HASk crystal structure at Hx = 0. However, energy of the states formed in the process of increasing Hx is apparently higher than that in the process of decreasing Hx, at the range of 0 < Hx < ~1.8. In this regard, we reckon that the alternating HSk and HASk crystal structure is a strongly frustrated spin system39. This is why the multi-step magnetization appearing in the <Mxy>-Hx curve in Fig. 6(a).


Before concluding this work, we briefly discuss experimental realization of the SC2 phase. B20-type helimagnets with another type of easy-plane anisotropy As(Siz)2 (As > 0) were also studied by MC simulation, which shows that this anisotropy may lead to the alternating HSk and HASk crystal structures on a square lattice at proper A too40. In addition, skyrmions in B20-type helimagnet thin films with easy-plane anisotropy As(Siz)2 (As > 0) were recently studied using a micromagnetic model including demagnetization and a three-dimensional geometry41. This work showed the demagnetization effects and/or the finite thickness effects on the skyrmion energetics and stability, which were explicitly demonstrated in earlier studies42,43. Therefore, B20-type helimagnet thin films seem to be appropriate for realizing the SC2 phase, because the uniaxial anisotropy (As), compass anisotropy (Ac) as well as the demagnetization can be altered significantly by varying the film thickness and the laterally confined geometries27,41,42,44. Furthermore, hedgehog skyrmions in chiral magnets with Rashba SOC were previously studied by variational energy calculations, which imposes a topological constraint for the center and boundary spins in a unit cell16. This constraint allows the hedgehog skyrmions but naturally exclude the HSk or HASk states in the chiral magnets. Therefore, whether the compass anisotropy can lead to the energetically stable HSk or HASk states in chiral magnets with Rashba SOC is yet to be explored.

On the other hand, temperature dependence of the Hall resistivity, susceptibility and magnetic modulation vector in MnGe was revealed by experiments11,13. These features are also important for our understanding of the underlying physics of MnGe, which would be an interesting challenge for our further studies.

In summary, we have presented in this work a detailed analysis of a classical microscopic spin model comprising FM exchange, DM interaction, compass anisotropy, and Zeeman energy, which may be widely used for the chiral magnets with strong SOC. We use classical MC simulation to determine the low temperature phase diagram as a function of compass anisotropy and out-of-plane magnetic field. In this phase diagram, we propose the alternating HSk and HASk crystal structures on a square lattice to be the candidate for the zero-field ground state of the helimagnet MnGe. The simulated evolution of the spin structure driven by magnetic field is in good accordance with experimental observations on MnGe. Therefore, our microscopic spin model successfully captures the main physics responsible for the magnetic structures in the helimagnet MnGe. This study may also provide a new theoretical guide to understand the magnetic states in other systems (such as the complex oxide heterostructure interfaces, and artificial ultracold Bose gas) with strong SOC, and may bring some insights for the future experiments.


Simulation algorithm

In this work, we perform MC simulations based on the standard single-flip Metropolis algorithm combined with the over-relaxation method. This clustering algorithm is believed to be effective in equilibrating the frustrated spin models with a large size45,46,47. Our unit MC step consists of one Metropolis sweep and ten over-relaxation sweeps. Most calculations are carried out on a 32 × 32 lattice with periodic boundary conditions, and then further checked on a lattice with bigger size of 64 × 64, in which case the simulation results remain the same. For reaching the ground state, an annealing simulation scheme is employed. First, a paramagnetic lattice is chosen as the initial lattice and each simulation cycle starts from sufficiently high temperature (T). Then the lattice is cooled down gradually until a very low T = 0.018,31. The as-obtained spin state is treated approximately as the ground state.

For tracking the evolution of non-zero field spin states, a ladder protocol following earlier reports39,48,49 is employed. The corresponding zero-field ground state is chosen as the initial lattice, and then magnetic field H is varied following the linear protocol H(j) = H(j −1) + δt, where H(j) (j = 1, 2, 3, …) is the field at stage j, δ is a constant, and t is time measured in MC steps. In the protocol, the lattice at H(j −1) is taken as the initial lattice for calculating the state at H(j). Besides, we discard the data in the first Ne MC steps needed for equilibration, and calculate the averages over the following Nc MC steps at every H(j) stage. Typically, both Ne and Nc are 5 × 105. During this procedure, field H is increased from H = 0 till a sufficiently big magnitude so that the lattice reaches the field-driven FM state. Subsequently, the field is reduced back to H = 0 for exploring potential hysteresis behavior if any.

Additional Information

How to cite this article: Chen, J. P. et al. Exotic skyrmion crystals in chiral magnets with compass anisotropy. Sci. Rep. 6, 29126; doi: 10.1038/srep29126 (2016).


We thank P. Chu for stimulating discussions. This work was supported by the National Key Research Programme of China (Grant No. 2016YFA0300101), the Natural Science Foundation of China (11234005 and 11374147), the Priority Academic Program Development of Jiangsu Higher Education Institutions, China, the Foundation for Fostering the Scientific and Technical Innovation of Guangzhou University, the Science Research Foundation for Young Teachers in South China Normal University (15KJ16), and the Foundation for Distinguished Young Talents in Higher Education of Guangdong (2015KQNCX023).


Author Contributions J.M.L. conceived the research project and J.P.C. performed the computations. J.P.C., D.W.Z., and J.M.L. commented the modeling and discussed the results. J.P.C. and J.M.L. wrote the paper.


  • Qi X. L. & Zhang S. C. The quantum spin Hall effect and topological insulators. Phys. Today 63, 33–38 (2010).
  • Hasan M. Z. & Kane C. L. Colloquium: topological insulators. Rev. Mod. Phys. 82, 3045 (2010).
  • Ishikawa Y., Tajima K., Bloch D. & Roth M. Helical spin structure in manganese silicide MnSi. Solid State Commun. 19, 525–528 (1976).
  • Mühlbauer S. et al. . Skyrmion lattice in a chiral magnet. Science 323, 915–919 (2009). [PubMed]
  • Bak P. & Jensen M. H. Theory of helical magnetic structures and phase transitions in MnSi and FeGe. J. Phys. C: Solid State Phys. 13, L881 (1980).
  • Yu X. Z. et al. . Near room-temperature formation of a skyrmion crystal in thin-films of the helimagnet FeGe. Nat. Mater. 10, 106–109 (2011). [PubMed]
  • Beille J., Voiron J. & Roth M. Long period helimagnetism in the cubic B20 FexCo1−xSi and CoxMn1−xSi alloys. Solid State Commun. 47, 399–402 (1983).
  • Yu X. Z. et al. . Real-space observation of a two-dimensional skyrmion crystal. Nature 465, 901–904 (2010). [PubMed]
  • Seki S., Yu X. Z., Ishiwata S. & Tokura Y. Observation of skyrmions in a multiferroic material. Science 336, 198–201 (2012). [PubMed]
  • Nagaosa N. & Tokura Y. Topological properties and dynamics of magnetic skyrmions. Nat. Nanotech. 8, 899–911 (2013). [PubMed]
  • Kanazawa N. et al. . Large topological Hall effect in a short-period helimagnet MnGe. Phys. Rev. Lett. 106, 156603 (2011). [PubMed]
  • Neubauer A. et al. . Topological Hall effect in the A phase of MnSi. Phys. Rev. Lett. 102, 186602 (2009). [PubMed]
  • Kanazawa N. et al. . Possible skyrmion-lattice ground state in the B20 chiral-lattice magnet MnGe as seen via small-angle neutron scattering. Phys. Rev. B 86, 134425 (2012).
  • Stishov S. M. & Petrova A. E. Itinerant helimagnet MnSi. Phys.-Usp. 54, 1117–1130 (2011).
  • Banerjee S., Erten O. & Randeria M. Ferromagnetic exchange, spin-orbit coupling and spiral magnetism at the LaAlO3/SrTiO3 interface. Nat. Phys. 9, 626–630 (2013).
  • Banerjee S., Rowland J., Erten O. & Randeria M. Enhanced stability of skyrmions in two-dimensional chiral magnets with Rashba spin-orbit coupling. Phys. Rev. X 4, 031045 (2014).
  • Cole W. S., Zhang S., Paramekanti A. & Trivedi N. Bose-Hubbard models with synthetic spin-orbit coupling: Mott insulators, spin textures, and superfluidity. Phys. Rev. Lett. 109, 085302 (2012). [PubMed]
  • Zhang D. W., Chen J. P., Shan C. J., Wang Z. D. & Zhu S. L. Superfluid and magnetic states of an ultracold Bose gas with synthetic three-dimensional spin-orbit coupling in an optical lattice. Phys. Rev. A 88, 013612 (2013).
  • Farrell A. & Pereg-Barnea T. Strong coupling expansion of the extended Hubbard model with spin-orbit coupling. Phys. Rev. B 89, 035112 (2014).
  • Nussinov Z. & van den Brink J. Compass models: Theory and physical motivations. Rev. Mod. Phys. 87, 1 (2015).
  • van den Brink J. Orbital-only models: ordering and excitations. New J. Phys. 6, 201 (2004).
  • Jackeli G. & Khaliullin G. Mott insulators in the strong spin-orbit coupling limit: from Heisenberg to a quantum compass and Kitaev models. Phys. Rev. Lett. 102, 017205 (2009). [PubMed]
  • Yildirim T., Harris A. B., Entin-Wohlman O. & Aharony A. Symmetry, spin-orbit interactions, and spin anisotropies. Phys. Rev. Lett. 73, 2919 (1994). [PubMed]
  • Bert J. A. et al. . Direct imaging of the coexistence of ferromagnetism and superconductivity at the LaAlO3/SrTiO3 interface. Nat. Phys. 7, 767–771 (2011).
  • Caviglia A. D. et al. . Tunable Rashba spin-orbit interaction at oxide interfaces. Phys. Rev. Lett. 104, 126803 (2010). [PubMed]
  • Butenko A. B., Leonov A. A., Rößler U. K. & Bogdanov A. N. Stabilization of skyrmion textures by uniaxial distortions in noncentrosymmetric cubic helimagnets. Phys. Rev. B 82, 052403 (2010).
  • Johnson M. T., Bloemen P. J. H., den Broeder F. J. A. & de Vries J. J. Magnetic anisotropy in metallic multilayers. Rep. Prog. Phys. 59, 1409 (1996).
  • King P. D. C. et al. . Large tunable Rashba spin splitting of a two-dimensional electron gas in Bi2Se3. Phys. Rev. Lett. 107, 096802 (2011). [PubMed]
  • Li X., Liu W. V. & Balents L. Spirals and skyrmions in two dimensional oxide heterostructures. Phys. Rev. Lett. 112, 067202 (2014). [PubMed]
  • Landau D. P. & Binder K. In A Guide to Monte Carlo Simulations in Statistical Physics 2nd edn, 68–190 (Cambridge, 2008).
  • Yi S. D., Onoda S., Nagaosa N. & Han J. H. Skyrmions and anomalous Hall effect in a Dzyaloshinskii-Moriya spiral magnet. Phys. Rev. B 80, 054416 (2009).
  • Makarova O. L. et al. . Neutron diffraction study of the chiral magnet MnGe. Phys. Rev. B 85, 205205 (2012).
  • Shibata K. et al. . Towards control of the size and helicity of skyrmions in helimagnetic alloys by spin-orbit coupling. Nat. Nanotech. 8, 723–728 (2013). [PubMed]
  • Binz B. & Vishwanath A. Chirality induced anomalous-Hall effect in helical spin crystals. Physica B 403, 1336–1340 (2008).
  • Park J.-H. & Han J. H. Zero-temperature phases for chiral magnets in three dimensions. Phys. Rev. B 83, 184406 (2011).
  • Lee S. Y. & Han J. H. Zero-energy bound states in a nodal topological lattice. Phys. Rev. B 91, 245121 (2015).
  • Alba E., Fernandez-Gonzalvo X., Mur-Petit J., Pachos J. K. & Garcia-Ripoll J. J. Seeing topological order in time-of-flight measurements. Phys. Rev. Lett. 107, 235301 (2011). [PubMed]
  • Pfleiderer C. Magnetic order: Surfaces get hairy. Nat. Phys. 7, 673–674 (2011).
  • Chen J. P. et al. . Manipulation of magnetic state in nanostructures by perpendicular anisotropy and magnetic field. J. Appl. Phys. 115, 243910 (2014).
  • Lin S. Z., Saxena A. & Batista C. D. Skyrmion fractionalization and merons in chiral magnets with easy-plane anisotropy. Phys. Rev. B 91, 224407 (2015).
  • Vousden M. et al. . Skyrmions in thin films with easy-plane magnetocrystalline anisotropy. Appl. Phys. Lett. 108, 132406 (2016).
  • Beg M. et al. . Ground state search, hysteretic behaviour, and reversal mechanism of skyrmionic textures in confined helimagnetic nanostructures. Sci. Rep. 5, 17137 (2015). [PMC free article] [PubMed]
  • Rybakov F. N., Borisov A. B. & Bogdanov A. N. Three-dimensional skyrmion states in thin films of cubic helimagnets. Phys. Rev. B 87, 094424 (2013).
  • Karhu E. A. et al. . Chiral modulations and reorientation effects in MnSi thin films. Phys. Rev. B 85, 094429 (2012).
  • Creutz M. Overrelaxation and Monte Carlo simulation. Phys. Rev. D 36, 515 (1987). [PubMed]
  • Campos I., Cotallo-Aban M., Martin-Mayor V., Perez-Gaviro S. & Tarancon A. Spin-glass transition of the three-dimensional Heisenberg spin glass. Phys. Rev. Lett. 97, 217204 (2006). [PubMed]
  • Viet D. X. & Kawamura H. Numerical evidence of spin-chirality decoupling in the three-dimensional Heisenberg spin glass model. Phys. Rev. Lett. 102, 027202 (2009). [PubMed]
  • Carubelli M. et al. . Spin reorientation transition and phase diagram of ultrathin ferromagnetic films. Phys. Rev. B 77, 134417 (2008).
  • Chen J. P. et al. . Stripe-vortex transitions in ultrathin magnetic nanostructures. J. Appl. Phys. 113, 054312 (2013).

Articles from Scientific Reports are provided here courtesy of Nature Publishing Group