---
author:
- 'I. Kézsmárki'
- 'S. Bordács'
- 'P. Milde'
- 'E. Neuber'
- 'L. M. Eng'
- 'J. S. White'
- 'H. M. Rønnow'
- 'C. D. Dewhurst'
- 'M. Mochizuki'
- 'K. Yanai'
- 'H. Nakamura'
- 'D. Ehlers'
- 'V. Tsurkan'
- 'A. Loidl'
title: |
    Néel-type Skyrmion Lattice with Confined Orientation\
    in the Polar Magnetic Semiconductor GaV$_4$S$_8$ 
---

**Following the early prediction of the skyrmion lattice (SkL) [@Bogdanov1989; @Bogdanov1989_2; @Bogdanov1994; @Bogdanov1994_2; @Rossler2006]—a periodic array of spin vortices—it has been observed recently in various magnetic crystals mostly with chiral structure [@Pfleiderer2004; @Muhlbauer2009; @Munzer2010; @Yu2010; @Yu2011; @Wilhelm2011; @Heinze2011; @Seki2012; @Adams2011; @Adams2012; @Tomomura2012; @Yu2012; @Romming2013; @Milde2013; @Shibata2013; @Yu2014; @Park2014; @Janson2014]. Although non-chiral but polar crystals with C$_{nv}$ symmetry were identified as ideal SkL hosts in pioneering theoretical studies [@Bogdanov1989; @Bogdanov1989_2; @Bogdanov1994; @Bogdanov1994_2], this archetype of SkL has remained experimentally unexplored. Here, we report the discovery of a SkL in the polar magnetic semiconductor GaV$_4$S$_8$ with rhombohedral (C$_{3v}$) symmetry and easy axis anisotropy. The SkL exists over an unusually broad temperature range compared with other bulk crystals [@Pfleiderer2004; @Muhlbauer2009; @Munzer2010; @Yu2011; @Seki2012; @Adams2012; @Milde2013] and the orientation of the vortices is not controlled by the external magnetic field but instead confined to the magnetic easy axis. Supporting theory attributes these unique features to a new non-chiral or Néel-type of SkL [@Bogdanov1989; @Bogdanov1994; @Leonov2014] describable as a superposition of spin cycloids in contrast to the Bloch-type SkL in chiral magnets described in terms of spin helices [@Pfleiderer2004; @Muhlbauer2009; @Munzer2010; @Yu2010; @Yu2011; @Wilhelm2011; @Seki2012; @Adams2012; @Tomomura2012; @Milde2013; @Shibata2013; @Park2014].**

In non-centrosymmetric crystals, the energy associated with ferromagnetic domain walls can be negative, hence, the homogeneous ferromagnetic state becomes unstable against SkL formation [@Bogdanov1989; @Bogdanov1989_2; @Bogdanov1994; @Bogdanov1994_2; @Rossler2006; @Yi2009]. Two basic types of magnetic domain walls can form different skyrmionic spin textures [@Bogdanov1989; @Bogdanov1994; @Bogdanov1994_2; @Leonov2014]. Bloch-type domain walls, where the spins rotate in the plane parallel to the domain boundary, can form whirlpool-like skyrmions with a given handedness as observed in chiral magnets [@Pfleiderer2004; @Muhlbauer2009; @Munzer2010; @Yu2010; @Yu2011; @Wilhelm2011; @Seki2012; @Adams2012; @Tomomura2012; @Milde2013; @Shibata2013; @Park2014]. In contrast, Néel-type domain walls, where the spins rotate in a plane perpendicular to the domain boundary, can produce non-chiral skyrmions. (See Fig. 5 for an illustration of the two skyrmion types.) This second skyrmion archetype, expected to emerge in polar magnets with C$_{nv}$ crystal symmetry [@Bogdanov1989; @Bogdanov1994; @Bogdanov1994_2; @Leonov2014], has not been observed yet. We reveal the formation of such a SkL in GaV$_4$S$_8$, a magnetic semiconductor from this crystal symmetry class, via magnetic susceptibility, atomic force microscopy (AFM) and small-angle neutron scattering (SANS) measurements.

GaV$_4$S$_8$, a member of the lacunar spinel family [@Phouc2013; @Elmeguid2004; @Dorolti2010; @Kim2014; @Guiot2013; @Singh2014; @Pocha2000], has a non-centrosymmetric cubic (T$_d$) structure at room temperature [@Pocha2000]. Lacunar describes the lack of every second Ga atom compared to the normal spinel structure. These ordered defects break the vanadium pyrochlore lattice into a network of alternating larger and smaller V$_4$ tetrahedra, the so-called breathing pyrochlore lattice. The magnetic building blocks are the smaller V$_4$ clusters with spin 1/2, and form a face centered cubic (FCC) lattice. The hybridization on a single V$_4$ unit leads to one unpaired electron occupying a triply degenerate cluster orbital. This orbital degeneracy is lifted by a cooperative Jahn-Teller distortion that drives a cubic to rhombohedral structural transition at $T_s$=42K [@Pocha2000]. In the polar rhombohedral (C$_{3v}$) phase the FCC lattice is stretched along one of the four cubic $\langle$111$\rangle$ axes [@Pocha2000].

![image](Figure1_final.pdf){width="6.8in"}

The structural transition creates a multi-domain state with submicron-thick sheets of the four different rhombohedral domains, as seen in Fig. 1a. At $T_s$ the magnetic exchange interaction changes from antiferro- to ferromagnetic and the material undergoes a magnetic transition at $T_C$=13K, a temperature very close to the Curie-Weiss temperature in the rhombohedral phase [@Yadav2008]. In the inset of Fig. 1b, the dependence of the magnetization curves on the field orientation indicate a considerable anisotropy. Indeed, our mean-field analysis shows the ground state to be an easy axis ferromagnet with $J_{\parallel}$$\approx$8.7K and $J_{\perp}$$\approx$8.3K as described in the Supplementary Information. Here, $J_{\parallel}$ and $J_{\perp}$ respectively denote the exchange coupling for spin components parallel and perpendicular to the easy axis, which itself coincides with the direction of the rhombohedral distortion. These exchange values are also consistent with $T_C$=13K.

GaV$_4$S$_8$ was considered as an ordinary ferromagnet below $T_C$ [@Pocha2000; @Yadav2008], though there is a report about low-field magnetic anomalies in this material [@Nakamura2009]. We found that the magnetic state is far from ordinary and contains several neighbouring phases. Below $T_C$, the low-field magnetization displays a sequence of steps indicating metamagnetic transitions (see Fig. 1b). The positions of the steps, associated with sharp peaks in the field-derivative of the magnetization, are different for magnetic field, $\mathbf{B}$, applied parallel to the \[111\], \[110\] and \[100\] cubic axes. The corresponding phase diagrams are shown in Figs. 1d-f. Similar low-field steps have been observed to separate the helical, SkL and conical states in MnSi [@Thessieu1997; @Lamago2006; @Muhlbauer2009], Mn$_{1-x}$(Co,Fe)$_{x}$Si [@Pfleiderer2010], FeGe [@Wilhelm2011] and Cu$_2$OSeO$_3$ [@Seki2012]. In GaV$_4$S$_8$ we assign the three different magnetic phases as cycloidal, SkL and ferromagnetic based on AFM imaging, SANS experiments and theoretical calculations as described below.

![image](Figure2_final.pdf){width="6.7in"}

For $\mathbf{B}$$\parallel$\[111\] there exist two additional magnetic states besides the cycloidal, SkL and ferromagnetic states. This arises due to multi-domain nature of the crystal and the easy axis anisotropy, since only one of the domains has its easy axis parallel to $\mathbf{B}$, while the easy axes of the other three domains—\[11$\bar{1}$\], \[1$\bar{1}$1\] and \[$\bar{1}$11\]—lie at 71 degrees to $\mathbf{B}$. For these latter three domains the non-collinear spin structures survive to higher fields because the magnetic anisotropy energy is much larger than the Zeeman energy, and only the field component parallel to the easy axis of each domain can influence its magnetic state. This leads to the two additional phases persisting up to $B$$=$80mT and 160mT, marked with asterisks in Fig. 1d, and assigned as the cycloidal and SkL states common for these three domains. For $\mathbf{B}$$\parallel$\[110\], the cycloidal and SkL phases common for the two domains with easy axes \[111\] and \[11$\bar{1}$\] are indicated in Fig. 1e. The other two domains with easy axes perpendicular to $\mathbf{B}$ cannot contribute to the magnetization is this field range. For $\mathbf{B}$$\parallel$\[100\], the easy axis of each domain span the same angle to $\mathbf{B}$, hence all domains share common phase boundaries as shown in Fig. 1f. Phase boundaries for all field orientations can be scaled together by projecting the field onto the easy axes of the corresponding domains (see Fig. 1c). Surprisingly, the SkL state exists over an unusually broad temperature range, down to $\sim$0.68$T_C$, compared with other bulk crystals where this phase is stable only over a 2-3% range immediately below $T_C$ [@Muhlbauer2009; @Seki2012; @Yu2011].

To observe directly the spin pattern in real space we performed AFM imaging in magnetic fields applied perpendicular to the sample surface using magnetic cantilever tips [@Milde2013]. For $\mathbf{B}$$\parallel$\[100\], Figs. 2a and b show images of the same surface area taken at $T$$=$8.5K and 8.9K in $B$$=$20mT and 50mT, respectively. The trail running along the \[001\] axis and dividing each image into two halves is a structural domain boundary.

![image](Figure3_final.pdf){width="7in"}

Although all four domains should host the same magnetic states for $\mathbf{B}$$\parallel$\[100\], we found that the magnetic pattern always changed at the domain boundaries, and often it was not possible to achieve good contrast for the two sides simultaneously. Since the wave vectors of the magnetic modulation ($\mathbf{q}$-vector) are expected to be perpendicular to the easy axis, they can be different for different domains, resulting in a rotation of the magnetic pattern between the two sides. In Fig. 2a, a modulated structure with a single $\mathbf{q}$-vector, representative of the zero- and low-field region, is observed with good contrast on the left side. Our calculations described later reproduce this phase as a cycloidal spin state, where the spins rotate in a plane containing the $\mathbf{q}$-vector. This is in contrast to the helical order generally found in other SkL materials, where the spins rotate in the plane perpendicular to $\mathbf{q}$. Assuming that the $\mathbf{q}$-vector is perpendicular to the easy axis, which is either the \[11$\bar{1}$\] or \[1$\bar{1}1$\] axis for the left-side domain, the periodicity of the cycloid in $B$$=$20mT is $a_{cyc}$$=$17.7$\pm$0.4nm. (For more details see the Supplementary Information.) This value is confirmed by our SANS study.

With increasing field the cycloidal pattern transforms into a distorted triangular lattice as shown in Fig. 2b. The same pattern rotated by about 90 degrees appears on both sides of the domain wall, though the contrast is weaker on the left side. As demonstrated in Fig. 2c, re-scaling the lattice on the right side by $\sqrt{3}$ perpendicular to the direction of the shortest lattice constant leads to a regular triangular lattice with a lattice constant of $a_{sky}$$=$22.2$\pm$1nm. This value is also consistent with our SANS data. The observed distortion is naturally explained by the orientational confinement of the SkL. Since the vortex cores are not aligned parallel to the magnetic field but instead the easy axis, which is either the \[111\] or \[$\bar{1}11$\] axis for the right-side domain, they form a regular triangular lattice in the hard plane. Therefore, the terminus of the vortices on the (100) surface form a distorted triangular lattice as sketched in Fig. 2d. When further increasing the field, the magnetic pattern disappears indicating a uniform ferromagnetic state.

SANS was employed to confirm the bulk nature of the magnetic structures and determine the $\mathbf{q}$-vectors directly. Sharp magnetic Bragg peaks were always observed at $\mathbf{q}$-vectors parallel to the cubic $\langle$110$\rangle$ directions irrespective of the field orientation, as shown in Fig. 3. For $\mathbf{B}$$\parallel$\[111\] the $\pm\mathbf{q}$-vectors form a regular hexagon, while for $\mathbf{B}$$\parallel$\[100\] a square of four spots is observed. These results, in agreement with the results of AFM imaging, confirm the orientational confinement of the SkL, i.e. the alignment of the vortex lines parallel to the easy axis of each domain. It also shows that in both the cycloidal and SkL states, the three possible $\pm\mathbf{q}$ pairs are parallel to the three \[110\] axes lying in the hard plane. For $\mathbf{B}$$\parallel$\[111\] a pale ring structure appears besides the hexagon in low fields. This implies that the $\mathbf{q}$-vectors are not fully pinned to the \[110\] axes but that they can point in any direction within the hard plane, at least in the surface and domain boundary regions observed by AFM.

![image](Figure4_final.pdf){width="6.4in"}

In each rhombohedral domain the cycloidal state can emerge with six different $\mathbf{q}$-vectors, while the SkL state is unique, and describable as a superposition of three cycloids whose $\mathbf{q}$-vectors sum to zero. First-order Bragg peaks observed using SANS do not allow a direct distinction between the cycloidal and the SkL states due to the presence of magnetic and structural domains. Nevertheless, the field dependence of the magnetic periodicity ($a=2\pi/|\mathbf{q}|$) clearly indicates the phase boundaries for both field orientations as shown in Fig. 3h. The coexistence of the cycloidal\* and SkL states in different domains for $\mathbf{B}$$\parallel$\[111\] is also traced via the splitting of $|\mathbf{q}|$.

![**$\mid$ Néel-type SkL in GaV$_4$S$_8$.** **a,** The Bloch-type skyrmion is a chiral object, i.e. it has no mirror-plane symmetry. **b,** In contrast, the Néel-type skyrmion has mirror-plane symmetry (combined with time reversal) for vertical planes through its center. **c,** Theoretically obtained temperature (T) versus magnetic field (B) phase diagram for the spin model by Eq. (\[eq:model\]) on the cubic lattice. The inset shows the direction of the Dzyaloshinskii-Moriya vectors for bonds within the $xy$ plane. The skyrmion lattice state corresponds to a Néel-type SkL.[]{data-label="Fig5"}](Figure5_final.pdf){width="3.2in"}

In GaV$_4$S$_8$, the S$=$1/2 spins sit on an FCC lattice stretched along the \[111\] direction, which can be regarded as triangular lattices stacked along \[111\]. To determine the spin patterns in the different phases we studied the following classical Heisenberg model on the triangular lattice by a Monte Carlo technique; $$\begin{aligned}
\mathcal{H}= \mathcal{H}_{\rm ex} + \mathcal{H}_{\rm DM} +
\mathcal{H}_{\rm Zeeman}, \label{eq:model}\end{aligned}$$ where $$\begin{aligned}
&&\mathcal{H}_{\rm ex}=-\sum_{<i,j>}\left(
 J_{\perp} m_{xi} m_{xj} + J_{\perp} m_{yi} m_{yj} +J_{\parallel} m_{zi} m_{zj} \right),\nonumber\\
%
&&\mathcal{H}_{\rm DM}
=\sum_{<i,j>}\bm D_{ij}(\bm m_i \times \bm m_{j}),\nonumber\\
%
&&\mathcal{H}_{\rm Zeeman}=-B \sum_{i}m_{zi}.\nonumber\end{aligned}$$ For the exchange interaction $\mathcal{H}_{\rm ex}$, we consider ferromagnetic coupling with XXZ-type anisotropy, where the $z$ (easy) axis is parallel to the \[111\] axis. The direction of the Dzyaloshinskii-Moriya vectors, $\bm D_{ij}$, are determined by the FCC crystal symmetry as shown in Fig. \[Fig4\]a. For the calculation we chose $J_{\parallel}/J_{\perp}$=1.06 and $D/J_{\perp}$=0.35, where $D=|\bm D_{ij}|$.

This model provides three ordered phases: the cycloidal, SkL and ferromagnetic states. Figure 4 displays the spin patterns obtained for the cycloidal and SkL states. The skyrmions exhibit the Néel-type domain wall alignment along the radial direction from their cores to their peripheries. Namely, the spins rotate within the plane parallel to the radial direction. The Fourier components of the magnetization configurations in the cycloidal and the SkL states are displayed in Figs. 4d and g, respectively. These results, in full agreement with SANS data, show that the $\pm\mathbf{q}$-vectors are parallel to the \[110\] axes in the hard plane in both states.

Extension of the model in Eq. (\[eq:model\]) to the distorted FCC lattice requires additional parameters as both $J_{\parallel}/J_{\perp}$ and $D/J_{\perp}$ can be different for bonds within and between the triangular lattice planes. Instead, we adapted the continuum model developed by Bogdanov and Hubert for $C_{nv}$ symmetry [@Bogdanov1994] to the cubic lattice in order to reproduce the main features of the phase diagram. Here, $\bm
D_{ij}$s are perpendicular to the bonds along the $x$ and $y$ directions as shown in the inset of Fig. 5c, and values of $J_{\parallel}/J_{\perp}$=1.06 and $D/J_{\perp}$=0.48 were used. We stress that irrespective of the details of possible models, GaV$_4$S$_8$ can only host a Néel-type SkL due to the orientation of Dzyaloshinskii-Moriya vectors governed by its C$_{3v}$ symmetry [@Bogdanov1989; @Bogdanov1989_2; @Bogdanov1994; @Bogdanov1994_2; @Leonov2014].

Figure 5c displays the phase diagram obtained for the simple cubic lattice. In agreement with experiment the triangular SkL appears only at finite temperatures. The skyrmions have the same Néel-type pattern, shown in Fig. 5b, as obtained for the triangular lattice model. The cycloidal phase remains stable even at $T$=0 in contrast to the experimental observation. This discrepancy likely arises from not including corrections due to quantum fluctuations in our classical Heisenberg model.

In a new class of materials we observed a SkL state emerging over an extraordinarily broad temperature range. While chiral or Bloch-type skyrmions have been investigated extensively in chiral magnets, lacunar spinels provide a unique arena to study another topological spin pattern, the non-chiral or Néel-type skyrmions. We found that the strong orientational confinement of the vortices ensures the robustness of two distinct skyrmionic states with a core magnetization pointing either up or down the easy axis. This may facilitate a unique magnetic control of the SkL, since we expect that the SkL induced by a magnetic field parallel to the easy axis can be rotated within the hard plane by an additional transverse field component. Such magnetic control is inconceivable in cubic helimagnets where vortex cores instantaneously co-align with the magnetic field. In addition, the polar crystal structure of lacunar spinels [@Singh2014; @Pocha2000] may be exploited for a non-dissipative electric field control of the SkL.

Methods {#methods .unnumbered}
=======

**Sample synthesis and characterization.** Single crystals of GaV$_4$S$_8$ were grown by the chemical vapour transport method using iodine as the transport agent. The crystals are typically cuboids or hexagonal slabs with masses ranging from 1-60mg. The sample quality was checked by powder X-ray diffraction, specific heat and magnetization measurements. The crystallographic orientation of the samples was determined by X-ray Laue and/or neutron diffraction prior to the magnetization, SANS and AFM studies. The magnetization measurements were performed using an MPMS from Quantum Design. The field dependence of the magnetization was measured in increasing temperature steps following an initial zero-field cooling to 2K. The temperature dependence of the magnetization was also measured in different fields. For both of these approaches, we found no hysteresis in the magnetization curves thus indicating that the SkL state cannot be quenched, i.e. stabilized as a metastable state, in this compound.

**Non-contact atomic force microscopy.** All AFM data presented were obtained by non-contact atomic force microscopy (nc-AFM) performed in an Omicron cryogenic ultra-high vacuum (UHV) STM/AFM instrument using the RHK R9-control electronics. The microscope features interferometric detection which as a side-effect allows the tuning of the cantilevers’ effective quality factor $Q_{eff}$ by varying the laser intensity used for the detection of the cantilever motion. For all measurements we used magnetically coated AFM-tips from the SSS-QMFMR series from NANOSENSORSin order to achieve sensitivity to magnetic forces.

In nc-AFM mode the frequency $f$ of an oscillating cantilever is measured via a phase locked loop and kept constant by the topography feedback loop. The frequency shift $\Delta f$ with respect to the resonance frequency $f_0$ far away from the sample provides a measure of the force gradient apparent between tip and sample. Additionally, the oscillation amplitude $A$ is kept constant by adjusting the excitation amplitude $A_{exc}$. Electrostatic interactions between tip and sample have been minimized using an additional FM-Kelvin control loop without any additional lateral information in the measured contact potential difference.

We could not obtain magnetic contrast in magnetic force microscopy mode, where the tip is retracted from the surface and the conservative magnetic forces are measured via the induced frequency shift $\Delta f$. Instead, a dissipative force interaction induced the magnetic contrast in the excitation channel when the tip was brought close to the sample surface. The nature of this dissipative interaction is subject to further investigations. We could obtain sufficient magnetic contrast only in images recorded below 10 K, which is likely related to the low spin density of the material corresponding to a single spin per formula unit and the strong fluctuations of the S=1/2 spins.

With knowledge of the cantilever’s spring constant $k$, the actual oscillation amplitude $A$, the resonance frequency $f$, the effective quality factor $Q_{eff,0}$ and the signal-to-drive ratio $R_0$ far away from the sample we compute the dissipated power from the measured excitation amplitude as $$P = \frac{\pi k f R_0A}{Q_{eff,0}}  \cdot A_{exc}.$$ For the images in figure 2 these parameters are $k=3.4$N$/$m, $A=3$nm, $Q_{eff,0}=4.28\cdot10^5$, $f=70.87$kHz, $R_0=30.92$$\mu$m$/$V with excitation amplitudes ranging from 150$\mu$V to 1.1mV.

**Small-angle neutron scattering.** Small-angle neutron scattering (SANS) was used to study the long-wavelength, microscopic magnetic states in a 25.5mg single crystal sample of GaV$_4$S$_8$. The SANS measurements were carried out using the D33 instrument at the Institut Laue-Langevin (ILL), Grenoble, France, and the SANS-II instrument at the Swiss Spallation Neutron Source (SINQ), Paul Scherrer Institut, Switzerland. In a typical instrument configuration, neutrons of wavelength 5[Å]{} were selected with a FWHM spread ($\Delta\lambda$/$\lambda$) of  10%, and collimated over a distance of 5.3m before the sample. The scattered neutrons were collected by a two-dimensional multi-detector placed 5m behind the sample.

The single crystal sample of GaV$_4$S$_8$ was mounted inside a horizontal field cryomagnet that was installed on the beamline so that the direction of applied magnetic field was approximately parallel to the neutron beam. The sample was oriented to have a horizontal scattering plane that included the \[100\] and \[111\] crystal directions. Applying the magnetic field along either of these directions was achieved by simply rotating the sample stick inside the cryomagnet. The SANS measurements were done by both tilting and rotating the sample and cryomagnet together through angular ranges that moved the magnetic diffraction peaks through the Ewald sphere. The SANS patterns presented in Fig. 3 were constructed by summing together the detector measurements taken at each rotation angle, which allows all of the Bragg spots to be presented in a single image.

**Theoretical calculations.** The Monte-Carlo calculations for Fig. \[Fig4\] are performed for the model in Eq. \[eq:model\] on the two-dimensional triangular lattice with $144  \times 144$ sites, while those for Fig. \[Fig5\] are for the model in Eq. \[eq:model\] on the three-dimensional cubic lattice with $72
\times 72 \times 8$ sites. For both cases periodic boundary conditions are imposed. The XXZ-type anisotropy in the exchange term comes from the orbital ordering with $3z^2-r^2$-type orbitals pointing in the \[111\] direction taken as the $z$ axis in the calculation [@Pocha2000]. The ratio of $J_{\parallel}/J_{\perp}$ was determined from a mean-field analysis of the low-temperature magnetization curves, with the value of $D/J_{\perp}$ fixed according to the observed cycloidal pitch, $a_{cyc}$. As long as the magnetic states of slow spatial variation such as cycloid and skyrmion are concerned, influence from complicated crystal structure becomes reduced, which justifies the choice to study the phase diagram using a lattice spin model on the simple cubic lattice. In the latter case, the Dzyaloshinskii-Moriya term is given by $$\begin{aligned}
\mathcal{H}_{\rm DM}=-D \sum_{i} \left[
 \hat{\bm y} \cdot (\bm m_i \times \bm m_{i+\hat{x}})
-\hat{\bm x} \cdot (\bm m_i \times \bm m_{i+\hat{y}}) \right].\end{aligned}$$ By comparing the energies of the three magnetic phases, we find that only the cycloidal and the ferromagnetic states are stable at $T$=0, while the SkL state is not. The transition points at finite $T$ are determined according to the calculated temperature-dependence of the specific heat and the net magnetization. The specific heat has a peak at the transition from the paramagnetic to the ordered phases. At the transition from the ferromagnetic phase to the cycloidal or the SkL phase, the magnetization shows a sharp drop with increasing temperature, while the specific heat shows almost no anomaly. Magnetic structures in these phases are identified from both the real-space and the reciprocal space mappings of their local magnetization configurations.

Bogdanov, A. N. & Yablonskii, D. A. Thermodynamically stable “vortices” in magnetically ordered crystals. The mixed state of magnets. [*Zh. Eksp. Teor. Fiz.*]{} [**95**]{}, 178-182 (1989).

Bogdanov, A. N. & Hubert, A. Thermodynamically stable magnetic vortex states in magnetic crystals. [*J. Magn. Magn. Mater.*]{} [**138**]{}, 255-269 (1994).

Bogdanov, A. N. & Hubert, A. The Properties of Isolated Magnetic Vortices. [*Phys. Stat. Sol. (b)*]{} [**186**]{}, 527-543 (1994).

Bogdanov, A. N. & Yablonskii, D. A. Contribution to the theory of inhomogeneous states of magnets in the region of magnetic-field-induced phase transitions. Mixed state of antiferromagnets. [*Zh. Eksp. Teor. Fiz.*]{} [**96**]{}, 253-260 (1989).

Rössler, U. K., Bogdanov, A. N. & Pfleiderer, C. Spontaneous skyrmion ground states in magnetic metals. [*Nature*]{} [**422**]{}, 797-801 (2006).

Pfleiderer, C. [*et al.*]{} Partial order in the non-Fermi-liquid phase of MnSi. Nature [**427**]{}, 227–-231 (2004).

Münzer, W. [*et al.*]{} Skyrmion lattice in the doped semiconductor Fe$_{1-x}$Co$_x$Si. [*Phys. Rev. B*]{} [**81**]{}, 041203(R) (2010).

Mühlbauer, S. [*et al.*]{} Skyrmion Lattice in a Chiral Magnet. [*Science*]{} [**323**]{}, 915-919 (2009).

Yu, X. Z. [*et al.*]{} Near room-temperature formation of a skyrmion crystal in thin-films of the helimagnet FeGe. [*Nat. Matter.*]{} [**10**]{}, 106-109 (2011).

Seki, S., Yu, X. Z., Ishiwata, S. & Tokura, Y. Observation of Skyrmions in a Multiferroic Material. [*Science*]{} [**336**]{}, 198-201 (2012).

Adams, T. [*et al.*]{} Long-Wavelength Helimagnetic Order and Skyrmion Lattice Phase in Cu$_2$OSeO$_3$. [*Phys. Rev. Lett.*]{} [**108**]{}, 237204 (2012).

Milde, P. [*et al.*]{} Unwinding of a Skyrmion Lattice by Magnetic Monopoles. [*Science*]{} [**340**]{}, 1076-1080 (2013).

Yu, X. Z. [*et al.*]{} Real-space observation of a two-dimensional skyrmion crystal. [*Nature*]{} [**465**]{}, 901–-904 (2010).

Tonomura, A. [*et al.*]{} Real-space observation of skyrmion lattice in helimagnet MnSi thin samples. [*Nano Lett.*]{} [**12**]{}, 1673-–1677 (2012).

Shibata, K. [*et al.*]{} Towards control of the size and helicity of skyrmions in helimagnetic alloys by spin–orbit coupling. [*Nat. Nanotechnol.*]{} [**8**]{}, 723–-728 (2013).

Park, H. S. [*et al.*]{} Observation of the magnetic flux and three dimensional structure of skyrmion lattices by electron holography. [*Nat. Nanotechnol.*]{} [**9**]{}, 337–-342 (2014).

Wilhelm, H. [*et al.*]{} Precursor Phenomena at the Magnetic Ordering of the Cubic Helimagnet FeGe. [*Phys. Rev. Lett.*]{} [**107**]{}, 127203 (2011).

Janson, O. [*et al.*]{} The quantum nature of skyrmions and half-skyrmions in Cu$_2$OSeO$_3$. [*Nat. Commun.*]{} [**5**]{}, 5376 (2014).

Heinze, S. [*et al.*]{} Spontaneous atomic-scale magnetic skyrmion lattice in two dimensions. [*Nat. Phys.*]{} [**7**]{}, 713-–718 (2011).

Adams, T. [*et al.*]{} Long-Range Crystalline Nature of the Skyrmion Lattice in MnSi. [*Phys. Rev. Lett.*]{} [**107**]{}, 217206 (2011).

Yu, X. Z. [*et al.*]{} Magnetic stripes and skyrmions with helicity reversals. [*PNAS*]{} [**109**]{}, 8856-8860 (2012).

Romming, N. [*et al.*]{} Writing and deleting single magnetic skyrmions. [*Science*]{} [**341**]{}, 636–-639 (2013).

Yu, X. Z. [*et al.*]{} Biskyrmion states and their current-driven motion in a layered manganite. [*Nat. Commun.*]{} [**5**]{}, 3198 (2014).

Leonov, A. Twisted, localized, and modulated states described in the phenomenological theory of chiral and nanoscale ferromagnets. Ph.D Thesis (2014).

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).

Ta Phouc, V. [*et al.*]{} Optical Conductivity Measurements of GaTa$_4$Se$_8$ Under High Pressure: Evidence of a Bandwidth-Controlled Insulator-to-Metal Mott Transition. [*Phys. Rev. Lett.*]{} [**110**]{}, 037401 (2013).

Abd-Elmeguid, M. M. [*et al.*]{} Transition from Mott Insulator to Superconductor in GaNb$_4$Se$_8$ and GaTa$_4$Se$_8$ under High Pressure. [*Phys. Rev. Lett.*]{} [**93**]{}, 126403 (2004).

Dorolti, E. [*et al.*]{} Half-Metallic Ferromagnetism and Large Negative Magnetoresistance in the New Lacunar Spinel GaTi$_3$VS$_8$. [*J. Am. Chem. Soc.*]{} [**132**]{}, 5704-5710 (2010).

Kim, H.-S., Im, J., Han, M. J. & Jin, H. Spin-orbital entangled molecular j$_{eff}$ states in lacunar spinel compounds. [*Nat. Commun.*]{} [**5**]{}, 3988 (2014).

Guiot, V. [*et al.*]{} Avalanche breakdown in GaTa$_4$Se$_{8-x}$Te$_x$ narrow-gap Mott insulators. [*Nat. Commun.*]{} [**4**]{}, 1722 (2013).

Singh, K. [*et al.*]{} Orbital-Ordering-Driven Multiferroicity and Magnetoelectric Coupling in GeV$_4$S$_8$. [*Phys. Rev. Lett.*]{} [**113**]{}, 137602 (2014).

Pocha, R., Johrendt, D. & Pöttgen, R. Electronic and Structural Instabilities in GaV$_4$S$_8$ and GaMo$_4$S$_8$. [*Chem. Mater.*]{} [**12**]{}, 2882-2887 (2000).

Yadav, C. S., Nigam, A. K. & Rastogi A. K. Thermodynamic properties of ferromagnetic Mott-insulator GaV$_4$S$_8$. [*Physica B*]{} [**403**]{}, 1474-1475 (2008).

Nakamura, H. [*et al.*]{} Low-field multi-step magnetization of GaV$_4$S$_8$ single crystal. [*J. Phys.: Conf. Ser.*]{} [**145**]{}, 012077 (2009).

Thessieu, C., Pfleiderer, C., Stepanov, A. N. & Flouquet, J. Field dependence of the magnetic quantum phase transition in MnSi. [*J. Phys.: Condens. Matter*]{} [**9**]{}, 6677-6687 (1997).

Lamago, D., Georgii, R., Pfleiderer, C. & Böni, P. Magnetic-field induced instability surrounding the A-phase of MnSi: Bulk and SANS measurements. [*Physica B*]{} [**385-386**]{}, 385-387 (2006).

Pfleiderer, C. [*et al.*]{} Skyrmion lattices in metallic and semiconducting B20 transition metal compounds. [*J. Phys.: Condens. Matter*]{} [**22**]{}, 164207 (2010).

**Acknowledgements** This work was supported by the Hungarian Research Funds OTKA K 108918, OTKA PD 111756 and Bolyai 00565/14/11, by the European Research Council Project CONQUEST, by the Swiss NSF Grant No. 200021\_153451, by the DFG under grant no. SFB 1143 and via the Transregional Research Collaboration TRR 80 From Electronic Correlations to Functionality (Augsburg/Munich/Stuttgart) and by JSPS KAKENHI under Grant Nos. 25870169 and 25287088 from MEXT Japan.

**Author Contributions** I.K., S.B., P.M., E.N., L.M.E., J.S.W., C.D.D., D.E., V.T. performed the measurements; I.K., S.B., P.M., E.N., H.M.R., J.S.W., A.L. analysed the data; V.T., H.N. contributed to the sample preparation; M.M., K.Y. developed the theory; I.K. wrote the manuscript and planned the project.

**Additional information** The authors declare no competing financial interests.
