Dzyaloshinskii-Moriya interactions, N\'eel skyrmions and V$_4$ magnetic clusters in multiferroic lacunar spinel GaV$_4$S$_8$

Using ab initio density functional theory with static mean-field correlations, we calculate the Heisenberg and Dzyaloshinskii-Moriya interactions (DMI) for an atomistic spin Hamiltonian for the lacunar spinel, GaV$_4$S$_8$. The parameters describing these interactions are used in atomistic spin dynamics and micromagnetic simulations. The magnetic properties of the lacunar spinel GaV$_4$S$_8$, a material well-known from experiment to host magnetic skyrmions of N\'eel character, are simulated with these ab initio calculated parameters. The Dzyaloshinskii-Moriya contribution to the micromagnetic energy is a sum of two Lifshitz invariants, supporting the formation of N\'eel skyrmions and its symmetry agrees with what is usually expected for $C_{3\nu}$-symmetric systems. The are several conclusions one may draw from this work. One concerns the quantum nature of the magnetism, where we show that the precise magnetic state of the V$_4$ cluster is crucial for understanding quantitatively the magnetic phase diagram. In particular we demonstrate that a distributed-moment state of each V$_4$ cluster explains well a variety of properties of GaV$_4$S$_8$, such as the band gap, observed Curie temperature and especially the stability of N\'eel skyrmions in the experimentally relevant temperature and magnetic-field range. In addition, we find that electronic correlations change visibly the calculated value of the DMI.


I. INTRODUCTION
Magnetic skyrmions, which are topological spin textures, have been found in a few materials classes (B20 compounds [1,2], Co-Mn-Zn alloys etc., see also a review in Ref. [3]) as well as in low-dimensional systems of transition metal multilayers (Pt/Co/Ta [4], Ir/Fe/Co/Pt [5], Pd/Fe/Ir(111) [6] etc.) where, for some systems, topological magnetism was observed even at room temperature and in the absence of applied magnetic field.Even more unique are bulk magnets where skyrmions coexist with ferroelectricity and the only known examples are Cu 2 OSeO 3 and the lacunar spinels GaV 4 S 8 [7], GaV 4 Se 8 [8,9], GaMo 4 S 8 [10] and GaMo 4 Se 8 [11].The spinel compounds are especially interesting, because the ferroelectricity is of rare orbital-driven origin and has a considerable magnitude and because the skyrmions are of Néel character.This is in contrast to all the other bulk systems, where only Bloch skyrmions are observed.
Such a unique behavior of lacunar spinels has been attributed to the C 3ν point group of the crystal structure (Fig. 1).It was argued [12] that this symmetry implies a specific form of the Lifshitz invariants describing the Dzyaloshinskii-Moriya interaction (DMI), which leads to the stability of Néel skyrmions and contrasts with the isotropic DMI in B20 compounds where only Bloch skyrmions emerge.At the same time, theoretical studies of the DMI in lacunar spinels are sparse.In many papers, some values of magnetic interactions between S = 1/2 V 4 -clusters are assumed and spin models with only nearest neighbors are used then to model the magnetic textures at varying external magnetic field strengths [7,9].On the other hand, there are a few works FIG. 1. a) Crystal structure of lacunar spinel GaV4S8 with V4 clusters as the main magnetic units with effective spin S = 1/2.The [111] axis (blue arrow) is shown, which corresponds to the C3 rotation axis of the C3ν point group.b,c) The V4 cluster of lacunar spinels in the low-temperature phase, where the distance between inequivalent V1 and V2 sites is 2.9 Å and d(V2 − V2) = 2.8 Å.The yellow arrows show the magnetic moments of individual V sites (m1 for V1 and m2 for V2) in b) the localized-moment ("L" state, m2 m1) and c) the distributed-moment ("D" state m2 ∼ m1) configurations.For details, see the main text.
where the Heisenberg and DM interactions between the V 4 -clusters are actually calculated, using perturbation theory or total-energy fitting methods [11,[13][14][15], and the results clarify the symmetry of DM vectors and the relative energy scales of different interactions in the system.Neutron experiments [16] and theoretical studies [17] propose that the magnetization is uniformly distributed over all V sites in the V 4 cluster.Nevertheless, it is not clear yet how the DMI in lacunar spinels is affected by electronic correlations and details of the magnetic state of the four-site transition metal clusters, which play the role of effective spins and form a face-centered network (Fig. 1).The work presented here is aimed at closing this gap by a systematic analysis of electronic and magnetic properties of the first skyrmionic lacunar spinel GaV 4 S 8 .

II. THEORETICAL METHODS
To understand the magnetic phenomena in GaV 4 S 8 we follow a multiscale approach where we start with a description on the level of individual electrons, proceed with atomistic magnetic interactions between effective V 4 cluster moments and, finally, model the magnetic textures at finite temperature in external field on length scales in the range 10 − 10 3 nm.Details about these three main steps are given further below.
Step 1 (electronic properties).The electronic structure and magnetic properties of lacunar spinels are studied in this work using density functional theory (DFT) [18] in the full-potential Linear Muffin-Tin Orbital implementation, available in the RSPt code [19,20].Electronic correlations are modeled here on the static mean-field level by means of the DFT+U approach with varied U and fixed Hund's coupling J H = 0.9 eV on top of the spin-polarized local-density or generalized-gradient approximations of the DFT.Summation in the Brillouin zone is performed on the shifted (16 × 16 × 16) k-mesh and the Fermi smearing with a temperature of 1 mRy is used for electronic occupations.All calculations are performed for the known experimental structure reported in the literature.
Two different V 4 cluster configurations are considered here, where the magnetic moment is either localized mostly on one V site or distributed over the whole cluster (four V sites, see Fig. 1c).Note that the total moment per cluster is around 1 µ B in both cases.It also has to be noted that, because of the elongation of V 4 tetrahedra along the [111]-direction, one of the V sites (V 1 ) is not symmetry-related to the other three sites (V 2 ) which are, however, equivalent to each other.As will become clear in the following, the cluster configuration can change the calculated properties substantially.
As our on-going calculations suggest and in accordance with literature [11], the other lacunar spinel GaMo 4 Se 8 with 4d states shows more uniformly magnetized Mo 4 clusters with parallel-aligned spin moments, in contrast to the 3d V-based spinels.This may indicate a fundamental difference between the 3d and 4d lacunar spinels, which will be discussed in a future work.
Magnetic interactions are calculated using the well-established Lichtenstein-Katsnelson-Antropov-Gubanov (LKAG) approach [21] (for a recent review see [22]), where the idea is to relate the interaction between two spins to the energy change due to a small perturbation of the magnetic state.We use the muffin-tin projection to calculate site-specific electronic parameters [23] and restrict ourselves to bilinear contributions to the magnetic energy for each pair of spins S i , which is written as follows: where one can distinguish between the isotropic Heisenberg (J ij ), Dzyaloshinskii-Moriya ( D ij ) and symmetric anisotropic ( Γij ) exchange interactions.To calculate the DM interaction, we perform three independent calculations where the spin axis is oriented along the x-, y-and z-directions to obtain the D x , D y and D z components.It is also worth mentioning that the total magnetic moment of GaV 4 S 8 varies only slightly upon such a global rotation of the magnetization.We also find that the symmetric anisotropic exchange Γij is one or two orders of magnitude smaller than the DM interaction for different bonds, so we do not include this type of exchange in further simulations for GaV 4 S 8 .Application of the whole approach described above to several transition metal systems in our previous works [24][25][26] has demonstrated the reliability of the calculated values of magnetic interactions, which justifies the use of this approach in the present work.
The critical temperature T c for the magnetic ordering is estimated from Monte Carlo simulations based on the calculated exchange interactions J ij and D ij .Simulations are done using the UppASD code [27,28] for bulk supercells containing (N ×N ×N ) unit cells with periodic boundary conditions, where we compared N = 10, 20, 30 to estimate the size effects (see example in Fig. 10a in the SI).Initial annealing is performed for 5 • 10 4 steps at the simulated temperature and statistical sampling of different observables is done afterwards for 10 5 steps.
To make it easier to discuss and present graphically different interactions in the studied spinels, we define effective magnetic interactions which characterize interactions between whole V 4 clusters (assuming frozen intracluster, magnetic degrees of freedom), instead of single atoms: Here, the summation runs over all sites of cluster m and all sites of another cluster n.Such effective parameters imply the assumption that the coupling between four spins within each V 4 cluster is significantly stronger than between different clusters, which is confirmed by our calculations, described below, and that during spin dynamics at not too high temperature the spins of the same cluster rotate synchronously.Numerical results obtained in this way are discussed below (Sections IV and V) for GaV 4 S 8 .These parameters, however, lead to a different temperature-dependent magnetization (from Monte-Carlo simulations) in the ferromagnetic state compared to the interactions J ij and D ij between individual V sites (see example in Fig. 10b in the SI).The nature of this effect will be studied in the future.
Step 3 (micromagnetics).From the atomistic interaction parameters J ij and D ij defining the spin model (1) one can go to the continuous limit and obtain the micromagnetic energy which can be used to model the magnetic properties on length scales that range up to hundreds of nanometers.Let us consider the derivation of the micromagnetic energy density ε DM due to the DM interaction (derivation for the Heisenberg exchange is similar).Following the derivations in [29,13], one can start from the atomistic DM interactions between spin on site i and all other spins on sites j: and replace S i with the micromagnetic order parameter, m ≡ m( r), as well as S j with 1 st -order expansion m + ( R ij • ∇) m, where R ij is the distance between the two sites.Substituting this into Eqn.(3) leads to a DM contribution to the micromagnetic energy density: where, one can define the spiralization matrix D ≡ D αβ (α, β = x, y, z) as In general, this matrix can contain nine non-zero components and has to be consistent with the crystal symmetry.For GaV 4 S 8 , we find that the only non-zero components are D xy = −D yx = D. Due to this specific form, the x-component of the vector in the square brackets in Eqn. ( 5) is D xy ∂/∂y and the y-component is −D xy ∂/∂x, while the z-component is zero.Accordingly, the DM contribution to the micromagnetic energy reads: A straightforward calculation gives the following energy density, ε DM , in a form which coincides with the interfacial type of DM interaction often discussed in the literature for magnetic films (see Eqn. (8) in Ref. 30): This result agrees with the Lifshitz invariants expected for the C 3ν crystal symmetry, as discussed, for example, in Ref. [12] (Eqn.6 in this reference) and Ref. [31] (Table I in this reference), and with the derivation for another spinel GaV 4 Se 8 in the SI of Ref. [13].In the latter reference, however, only nearest neighbors and cluster-cluster interactions were taken into account, while in our work we include the full information on the intersite interaction parameters J ij and D ij between several hundred and a couple of thousand neighbors.As a side remark, for bulk systems with cubic crystal symmetry and isotropic DMI, such as B20 compounds MnSi and FeGe, the second line in the determinant in Eqn.(7) would be D • (∂/∂x, ∂/∂y, ∂/∂z), since the spiralization matrix D is diagonal, as verified by our direct calculations (not shown here).This would lead to the usual expression for the isotropic DM energy density ε DM = D m • ( ∇ × m).Calculation of the micromagnetic parameters and magnetic textures (helical states and skyrmions) for B20 systems is discussed, for example, in our recent works [24,26] as well as in Refs.[32][33][34].
Similarly, one can also derive the magnetic energy for the Heisenberg exchange: Here, the first term is a constant energy contribution, since it is proportional to m 2 , that is equal to 1 (because constant length of magnetic moment vectors is considered).One can also show that the next term with a 1 st -order derivative reads: which is zero, because m α m α = 1.Finally, the leadingorder term in Eqn.(10) reads: Usually, only the diagonal term is considered in micromagnetic simulations (α = β), and the energy becomes proportional to the spin stiffness A defined as follows: As discussed in the literature [35], and our recent work [26], the numerical evaluation of micromagnetic parameters according to Eqns. ( 6) and (14) shows a convergence problem with respect to the real-space cutoff.For that reason, following the literature recipe in Ref. [35], we introduce an exponential regularization factor (in contrast to Eqn. (2), indices i and j refer either to atomic sites or different clusters): (15) where the limit µ → 0 is taken at the final step.The exponential factors are introduced here to improve the convergence with respect to the real-space cutoff for the magnetic interactions.The calculated values of A and D as functions of parameter 1.0 < µ < 2.0 are then extrapolated to µ = 0 using a 3 rd -order polynomial, which gives a reasonable fitting quality, and a discussion of some further technical details can be found in the SI of Ref. 26.
As mentioned on page 2, we compare the properties of GaV 4 S 8 for two electronic configurations of V 4 clusters which we find to be relatively close in energy (∆E ∼ 100 meV).In case of the distributed-moment configuration (Fig. 1b), we calculate the spin stiffness A and DM spiralization D from the effective interactions defined by Eqn. ( 2), because the spin stiffness from the original J ij interactions between individual V sites is negative.This is natural to expect since the cluster configuration is ferrimagnetic in our calculations.The effective cluster-cluster DM interaction, on the other hand, disregards the internal DMI between V sites in the same cluster, which should not matter for large-scale magnetic textures.This is because we assume that the internal magnetic structure of V 4 clusters is frozen, which is a good approximation due to the large calculated intracluster magnetic exchange which can have a magnitude as large as 28 meV.For the localized-moment state (Fig. 1c), we compute the micromagnetic parameters calculated from the original interatomic interactions while taking into account only the V 1 sites with the largest magnetic moment.The resulting micromagnetic parameters are scaled down by the square of the V 1 moment to simulate an effective 1 µ B V 4 cluster.
Regarding the DM spiralization matrix (Eqn.15), initially we obtain it in the basis of Cartesian unit vectors e α (xyz-basis).It turns out, however, that it is more convenient to work in the basis where the e z is along the [111] direction (relevant for skyrmions) and the other two basis vectors are in the (111)-plane.In particular, we choose e y as a vector connecting the [111] line with the end of the lattice vector a 2 and e x is obtained as a normalized cross product of e y and e z .The inversion of the matrix, where rows are formed by these new unit vectors e α , results in a unitary matrix Û that defines a transformation to the [111]-basis.For the DM spiralization matrix D in this new basis, one gets: where D is the matrix in the xyz-basis.
It is in the [111]-basis where the spiralization matrix of lacunar spinel GaV 4 S 8 has a dominant component D xy = −D yx = D, in agreement with the C 3ν crystal symmetry.For that reason, in the paper we discuss the DM interaction and perform the micromagnetic simulations in the basis where the z axis is along the [111]direction.
Regarding the on-site anisotropy, we assume in our simulations that the uniaxial anisotropy energy constant K 1 is between (10 − 16) kJ/m 3 , as suggested in previous works [36,37].This anisotropy constant is even smaller than K 1 = 45 kJ/m 3 for hcp Co and cubic anisotropy constant K 1 = 48 kJ/m 3 for bcc Fe but larger than K 1 = −0.5 kJ/m 3 for fcc Ni.Despite the small value of K 1 for GaV 4 S 8 , as we will see in the following, it is important to include anisotropy in simulations of magnetic skyrmions in this class of systems.
Using the calculated micromagnetic parameters, including the DM interaction (Eqn.8) and the literature values of the uniaxial anisotropy (10 − 16) kJ/m 3 , we perform micromagnetic simulations of magnetic textures at finite temperature and in external magnetic field using the UppASD [27,28] and MuMax3 [38] codes.The temperature is varied between 0 and 30 K and the external field -between 0 and 600 mT, which corresponds to the experimentally studied range of parameters.The micromagnetic region is described by a (512 × 512 × 1) mesh with an equidistant step ∆h in all directions (∆h = 0.5 nm for simulating skyrmions and ∆h = 1.0 nm for simulating ferromagnetic state and spin spirals) and nonperiodic boundary conditions.We have also verified the effect of dimension in the z-direction and the boundary conditions (see discussion in Section VI).The magnetization dynamics is described by the Landau-Lifshitz-Gilbert equation [39,40]: where m i is the magnetization of a given micromagnetic region i, and the effective field B i is determined from the micromagnetic parameters A, D and uniaxial anisotropy, K, as well as the external field, B, and the classical dipolar field; γ is the gyromagnetic ratio.At finite temperature, random field proportional to √ α T is added to B i , and the damping constant α is set to 0.04.Variable time step in the range 3 • 10 −15 − 5 • 10 −14 s is used for the dynamics simulations, and the total simulation time was around (2 − 5) ns, which allowed to reach the equilibrium state starting from a random magnetic configuration.
Micromagnetic simulations are run starting from a random magnetic configuration at a temperature of 20 K.The system is cooled in 2 K-temperature steps down to 0 K, which corresponds to simulated annealing, and then the system is heated up to 20 K with the same speed.Results for the anisotropy values K = 10 kJ/m 3 and K = 16 kJ/m 3 are compared.Since the [111]-axis of the crystal structure is defined now as the z-axis, which leads to expression (8) for the DMI energy, the anisotropy easy-axis is also along the z-direction.

Data Availability
The datasets used and/or analysed during the current study are available from the corresponding author on reasonable request.

III. ELECTRONIC PROPERTIES
In the following, we discuss the electronic structure, and related properties, of the lacunar spinel GaV 4 S 8 , calculated using the local-density (LDA) as well as the generalized-gradient (GGA) approximation, also including correlations described on the static mean-field level of spin-polarized DFT with Hubbard-U corrections.It is noteworthy that the electronic structure work presented here shows two stable minima, one with a total moment of ∼ 1 µ B for each tetrahedron of V atoms of the lacunar spinel structure (Fig. 1c) and one solution in which the total moment per V tetrahedron still is close to 1 µ B , but where one of the V atoms of each tetrahedron carries a significantly larger moment compared to the other V atoms (Fig. 1b).We refer to the first configuration as the distributed moment case, and the second configuration as the localized moment case.
The calculated electronic states of GaV 4 S 8 near the Fermi level, E F , are represented by a number of states with a narrow bandwidth crossing E F , as indicated by the LDA results for the localized moment state of V 4 clusters in Fig. 2a.The spin-splitting of bands increases somewhat in response to electronic correlations in DFT+U calculations and a small gap of the order of 0.1 eV appears in the electronic spectrum for U = 2 eV (Fig. 2c).Also, the distributed moment state (Fig. 1b) of V 4 clusters becomes stable, but we find that the energy of the localized moment state is somewhat lower by around 76 meV/f.u.Varying U between 0 and 2 eV increases the total magnetic moment per formula unit from 0.8 µ B to almost 1.0 µ B , and these values agree with the range of measured values in the literature [41].
The results discussed above are obtained within the local-density approximation.Using the generalizedgradient approximation (GGA) we obtain essentially the same trends, but GGA shows a stronger tendency to magnetism.This leads to the magnetic moment of each V 4 tetrahedron being close to 1.0 µ B already for pure DFT (U = 0 eV) calculations (Table I).At increasing correlation strength U , the moments of both V 1 and V 2 sites are enhanced (Table I), but the total moment of each V 4 cluster remains almost the same.Next, we notice that the DFT+U band structures and densities of states (Fig. 2) within LDA and GGA are similar but there is an offset around 1 eV in terms of the U values between the two approximations, meaning that in GGA+U one needs U values smaller by roughly 1.0 eV to get results similar to LDA+U .At U = 2 eV, we find again that the distributed-moment state can be stabilized (last row in Table I) but it is higher in energy by 99 meV/f.u.compared to the localized-moment state.Finally, the band gap calculated in GGA is larger compared to the LDA estimate, which is in accordance with literature (see Fig. 10 FIG. 2. Density of states of GaV4S8 calculated within density functional theory in LDA and GGA approximations without and with Hubbard-U corrections.Red and blue curves correspond to spin-up and -down states.The system is in the localized-moment semi-metallic state in GGA at U = 0 eV (no bands cross EF, but the energy gap is zero valued) and is semiconducting at U = 2 eV with a gap Eg = 0.2 eV (localized-moment state, dark lines) and 0.3 eV (distributedmoment state, light lines).
in [42]).TABLE I. Calculated magnetic moments of V1 and V2 sites and total moment for GaV4S8 with experimental structure.Results are obtained within GGA+U for the localizedmoment ("L" state) and distributed-moment ("D" state) configurations (Fig. 1b,c).At U = 0 eV, both in LDA and in GGA, we find a considerable ferromagnetic interlayer Heisenberg interaction (J ⊥ ) of the order of 4 K, which couples the neighboring (111)-planes of V 4 clusters.Furthermore, we find a much weaker intralayer interaction (J || ), which is between the neighboring clusters in each of these (111)planes.This picture changes when correlations are included within DFT+U .In particular, the intralayer ex- change J || becomes stronger (reaching up to almost 12 K) and can outweigh the interlayer exchange J ⊥ .Interestingly, within the LDA, J ⊥ increases as a function of U until U = 1 eV and then decreases, while GGA results show decreasing J ⊥ which even becomes antiferromagnetic at U = 2 eV where the localized-moment state of the V 4 clusters changes to the distributed-moment state (Fig. 1c).In LDA, however, such a transformation of the V 4 cluster does not lead to change of sign of effective magnetic interactions.
It is necessary to emphasize that the parameters J ⊥ and J || that we discuss here are actually the effective interactions defined by Eqn.(2) and are introduced to make the discussion of the results more transparent.For Monte Carlo (MC) simulations at different temperatures, we use the original data, i.e. interactions between individual V sites (even weakly magnetic ones).One should note that, even though the effective interlayer exchange becomes AFM for U = 2 eV in GGA, the Monte Carlo simulations using not just effective but all calculated intersite interactions actually find the correct ferromagnetic ground state, allowing to conclude that the effective cluster-cluster models probably cannot cover the whole physics in this system (see further Monte Carlo simulations in Fig. 10b in the SI).
Figure 3 shows how the Curie temperature, estimated from our MC simulations using both the calculated Heisenberg and Dzyaloshinskii-Moriya interactions, changes as a function of electronic correlation strength characterized by the U parameter.For the localized-moment state, the ordering temperature for U = (1.0 − 1.5) eV is around 12 K and agrees nicely with the measured value of 13 K [7].Electronic correlations increase slightly the Curie temperature (by a few degrees) and at some point (U = 2 eV) allow to sta- bilize the distributed moment state, as discussed above.The temperature dependence of magnetization (M (T )) curves for both cluster configurations is similar (Fig. 4) and look like typical M (T ) curves for ferromagnets.However, the distributed-moment state is characterized by a higher Curie temperature (around 23 K), which is overestimated compared to the experimental value of 13 K. Stronger magnetism for the distributed-moment configuration is likely caused by the larger number of magnetic exchange paths when all V sites in each V 4 cluster have non-zero magnetic moments.

V. DZYALOSHINSKII-MORIYA INTERACTION
For GaV 4 S 8 within pure DFT (LDA approach, U = 0 eV), we obtain the spiralization matrix in the [111]basis (in units of meV• Å) 43 , as described in the Methods section: The [111]-basis appears to be more convenient for studying the micromagnetic behavior than the Cartesian basis, since the D matrix has just one dominating component D xy which is to be substituted as the D parameter in Eqn.(8).The obtained form of the spiralization matrix agrees with a previous work [13] and the observation of Néel skyrmions in this bulk system [7].Possible origin of the slight asymmetry (∼ 10 −3 ) of the calculated spiralization matrix may be related to a weak dependence of electronic properties on the total magnetization direction due to spin-orbit coupling.
When electronic correlations are included on the meanfield level with LDA and U = 1 eV, the DMI increases dramatically by more than a factor of 4 and changes sign: Similar response to moderate correlations is observed for the cluster-cluster DM interactions, in particular, for the nearest-neighbor in-plane interaction D 1a (Fig. 5).Large enhancement of the DMI can be related to the change of the magnetic state of the V 4 clusters: at U = 1 eV the magnetic moment of one of V sites increases by more than a factor of 2, while the other three sites remain weakly magnetic but change to the opposite direction.This more asymmetric distribution of magnetization at U = 1 eV can, in principle, create additional inversion-symmetry breaking effect which increases the DMI.On top of that, the system becomes semimetallic (Fig. 2), since the band gap is zero but no bands directly cross the Fermi level.Within the generalized-gradient approximation (GGA), we find a contrasting behavior, since the DMI parameter D xy decreases in absolute value for the localized-moment state as a function of electronic correlations U , from D xy = −0.20 meV • Å at U = 0 eV to −0.16 meV • Å at U = 2 eV.We should mention that, in contrast to LDA, correlations added within GGA lead to a band gap opening already at U = 1 eV (Fig. 2) which could partially explain the differences in the calculated DMI.Notably, for the distributed-moment state in GGA, which is also stable at U = 2 eV, the DM parameter D xy = −1.37 meV • Å is considerably larger.Strong dependence of the DM spiralization constant on the cluster state makes sense in view of the other findings shown in Fig. 5 (for LDA approximation) where the atomistic DMI is larger for the more asymmetric magnetization distribution in the cluster which may break further the inversion symmetry.Also, the localized-and distributed-moment states of V 4 cluster lead to distinct band structures (Fig. 8 in the SI) and this can have a sizeable effect on the DMI too.
In addition, the spin stiffness A decreases by almost a factor of 5 (in GGA) in response to additional electronic correlations with U = 2 eV.Because of that, the A/D xy ratio decreases from 26.3 nm to 7.4 nm, suggesting again that the magnetic properties of GaV 4 S 8 are very sensitive to electronic correlations.On the other hand, at U = 2 eV also the distributed moment state is stable and it shows a relatively large spin stiffness, compared to the localized moment state.Also, the DM interaction is considerably larger resulting in the A/D ratio of 5.0 nm.The A/D ratio is important in the context of non-collinear magnetism as well as skyrmions and lower values of A/D, in general, are expected to indicate more compact skyrmions and helical spin states, which is confirmed by the micromagnetic simulations in Sec.VI.Notably, the LDA+U results show a spin stiffness which is roughly a factor of two larger compared to the GGA+U estimates for U ≥ 1 eV.
It is worth mentioning that the strongest interatomic DM interaction (D 0 ∼ 0.6 meV ≈ 6.8 K) in these spinels, according to our calculations, comes from the V-V bonds within each metal cluster, implying that the actual magnetic state of V 4 clusters may be non-collinear.On the other hand, the non-collinearity is not expected to be large, since the canting angles within the cluster should be of the order of D 0 /J 0 where J 0 is the nearest-neighbor Heisenberg V-V exchange which is, as we find, antiferromagnetic and in the range of several hundred Kelvin.For that reason, canting angles around a few degrees can be expected, which should not change the main findings reported in the present work.This intracluster DMI would actually contribute significantly to the calculated micromagnetic constant D, and it is a subtle question whether to include this DMI or not when addressing the behavior of skyrmions in this system.A strong argument not to do so, we suggest, is that the internal magnetic exchange in each V 4 cluster is very large in our calculations and, for that reason, the spins of each cluster are expected to co-rotate.In that case, the internal DMI as well as the internal Heisenberg exchange do not contribute to the micromagnetic energy.Table II summarizes our findings for the spin stiffness and DM spiralization for the different calculation setups and two V 4 cluster states.Å2 and spiralization Dxy in units of meV• Å) for GaV4S8 with experimental structure.The A and Dxy parameters are obtained from V1-V1 (in case of localizedmoment state, "L") or cluster-cluster (in case of distributedmoment state, "D") interactions, which exclude the intracluster exchange.A and Dxy for the "L" state are divided by the square of the magnetic moment of the V1 site.Estimate of the Curie temperature (TC) is obtained from classical Monte Carlo simulations using V-V interaction parameters for all V sites.Results obtained within GGA+U for different correlation strengths U and cluster states "L" and "D" (see main text and Fig. 1b,c Part of the motivation for this study is the experimental realization of Néel skyrmions in the lacunar spinel GaV 4 S 8 .As an ultimate test of the accuracy of the calculated electronic structure and interatomic exchange parameters, as well as the transition to micromagnetic interaction strengths, we explore here the possibility of skyrmion formation with the aforementioned parameters.It should be noted that so far no fitting to experimental data has been made and all calculations are made in abinitio mode.To undertake this investigation we compare three sets of calculated micromagnetic parameters and their ability to reproduce skyrmions.The parameters used are: Note that the values are different here when compared to those in Table II, because they are divided by the unit cell volume and converted to the units usually used in experimental reports.The saturation magnetization of the simulations is 41.35 kA/m which corresponds to 1 µ B per formula unit, a value found in experiment as well as in the calculations (Sec.III).
For the parameter set "c" (distributed-moment state), we obtain a multi-domain ferromagnetic state with isolated skyrmions (Fig. 6a) at zero temperature and zero  II).The color code shows the out-of-plane direction of the magnetic moment vector in each point of space.
applied field.The emergence of a spin-spiral magnetic state (Fig. 6b), for zero-field simulations at low but finite temperature, with a period of around 20 nm, agrees well with the measured value a cyc = 17.7 nm (see Fig. 3 in Ref. 7).In addition, a state with stable skyrmions, with calculated topological charge close to ±1, is found when the external magnetic field is applied along the ±zdirection with a strength between (50 − 300) mT (see Fig. 6c).The skyrmion size in our simulations depends on the external field and ranges from 27 nm at B = 25 mT to 13 nm at B = 300 mT, where the number of skyrmions is dramatically decreased.This estimated size is compatible with the experimentally observed skyrmion size a sky = 22.2 nm reported in Fig. 3 of Ref. 7. At higher fields, the system is ferromagnetic up to temperatures around (12 − 14) K.The latter marks the critical temperature also for other types of magnetic order in this system.All these findings are summarized in the calculated phase diagram in Fig. 7a and are in a good agreement with experiments in Ref. 44.
In contrast, the parameter sets "a" and "b" (obtained from the localized moment state) produce practically no skyrmions and show a magnetic order only at relatively low temperatures (see phase diagram in Fig. 7b), which is due to the smaller values of the A and D parameters.In general, larger spin stiffness for the distributed moment state can be explained by the fact that there are more interaction paths in the structure, since all V atoms are then magnetic.By varying the micromagnetic parameters A, D and K, we find (data not shown) that the A/D, A/K and D/K ratios are all important for stabilizing skyrmions, which explains why the parameter set "c" (distributed moment state) gives a better agreement with experiment, given that the anisotropy is in the range (10 − 16) kJ/m 3 .
Our conclusions for the three parameter sets used in the micromagnetic simulations are qualitatively robust with respect to at least 10%-variations of the strength of the interactions.We note that lower DM values lead to smaller saturation fields for stabilizing a ferromagnetic state.Results for non-periodic and periodic boundary conditions in the simulations are also found here to yield similar results.We have also verified that increasing the out-of-plane dimension (parallel to the external magnetic field; the z-axis) of the simulation cell from n z = 1 to n z = 16 does not change qualitatively the phase diagrams shown in Fig. 7 (which are obtained using n z = 1).The most prominent, quantitative change of these simulations, compared to the thin 2D simulation cell, is an increase of the Curie temperature for the distributed moment configuration up to around 20 K and a lowering of the magnetic field needed to induce the ferromagnetic state.For the localized-moment configuration, the Curie temperature remains essentially the same for thinner and thicker simulation cells.
To summarize this section, our results indicate that the distributed-moment configuration (Fig. 1c) describes better the magnetic phase diagram (Fig. 7a) of bulk GaV 4 S 8 spinel, which shows ferromagnetic phase, spin spirals and Néel skyrmions.An example of the magnetic texture of a Néel skyrmion is shown in Fig. 12, isolated in a ferromagnetic background a) and in a lattice b).Based on the experimental phase diagram reported in Ref. 44, the results shown in Fig. 7a) have similar trends for the transition to the paramagnetic phase.In addition, the transition to the ferromagnetic phase occurs at small magnetic fields (in the range around 250 − 400 mT-as shown in Fig. 7a)) which agrees with observations (in the experiments it is 50 − 160 mT).The Curie temperature calculated in our simulations is in the neighborhood of 17 K which is close to the experimental one, i.e. ∼ 13 K.The only difference between phase diagrams is related to the skyrmionic and cycloidal regions.In the experimental phase diagram, the skyrmion lattice phase only appears from 9 K to 12.7 K but in the theoretical phase diagram, the temperature range where both phases appear goes from 0 K until near 17 K.Figures 7a) and b) show the phase diagram of the system in 2D for distributed and localized moments, respectively.By performing simulations for the 3D system, we obtained the same Curie temperature for the localized moments as the 2D case but the phase transition to the ferromagnetic ordering was occurring at very small magnetic fields (around 15 − 20 mT).In the 3D simulations for the distributed moments, the calculations indicate that the same phase transitions for skyrmionic and cycloidal regions occur at the same values for the temperature and external magnetic field but the transition to the paramagnetic ordering arises at higher Curie temperature.In Fig. 13 in the SI are also shown eight layers of the 3D simulation for the distributed moments.The figure emphasizes the tubular quasi-2D nature of skyrmions in GaV 4 S 8 .

VII. CONCLUSIONS
In this theoretical study of the lacunar spinel GaV 4 S 8 , we find that the Dzyaloshinskii-Moriya interaction (DMI) calculated from first principles and the associated micro-magnetic energy reflect the C 3ν crystal symmetry and support the formation of Néel skyrmions, previously reported for this system [7].In contrast to previous works [13,14], we obtain a detailed picture of the magnetic interactions, both between individual V sites and between different V 4 clusters, not just the nearest neighbors.
Electronic correlations are important in this multiferroic system, since, on the one hand, they open up a semiconducting band gap ∼ 100 meV and, on the other hand, they enhance significantly the energy scale of the Dzyaloshinskii-Moriya interactions (spiralization D in Eqn. ( 6)) relative to the Heisenberg exchange (spin stiffness A in Eqn. ( 14)).In particular, within the generalized-gradient approximation we find a smaller A/D ratio of 7.4 nm for U = 2 eV (moderate correlations, localized-moment state) compared to the pure DFT result A/D = 26.3nm, even though the absolute value of D is reduced by correlations.We believe that this behavior is related to the opening of the band gap and redistribution of the magnetic density in the V 4 cluster.Compared to the localized-moment state, the ratio A/D = 5.0 nm for the distributed-moment configuration is remarkably smaller, while the spin stiffness A is a factor of 6 larger, allowing to distinguish these two cluster states based on the predicted properties.
Based on our micromagnetic simulations using our computed first-principles parameters, we conclude that a small |A/D| ratio is important for stabilizing Néel skyrmions with a size 13 − 27 nm close to the measured one (a sky = 22.2 nm, [7]), while the value of the spin stiffness A determines the critical temperature for magnetic order.Although the uniaxial magnetic anisotropy is weak, it has a considerable effect on the magnetic phase diagram (Fig. 7).With the literature values of the anisotropy, we find that the distributed-moment configuration, where all four V sites in each V 4 cluster have sizeable moments, describes much better the magnetic properties and textures in GaV 4 S 8 for the experimentally studied temperature and magnetic field range.
We note that there is a difference between the atomistic and micromagnetic results for this system.For example, in Fig. 4 the Curie temperature, T C , of the two different types of electronic (and magnetic) configurations are 15 K and 23 K, while the respective estimates from Fig. 7 are 4 K and 15 K, meaning a shift around 10 K between the 3D atomistic and 2D micromagnetic results.We have also made 3D micromagnetic simulations, which showed somewhat larger T C for the distributed-moment cluster, but the same T C for the localized-moment cluster.Overall, it is from these types of simulations difficult to pinpoint an ordering temperature with an accuracy of a few kelvin, 28 so that from these comparisons it is difficult to identify which electronic configuration (distributed or localized moments) is relevant for this system.However, magnetic textures at lower temperatures are more faithfully reproduced by the type of calculations/simulations presented here 28 , and for this reason we argue that the four-site, distributed moment configuration in GaV 4 S 8 , which is crucial to reproduce the magnetic properties, represents the correct electronic configuration of this material.
It should be noted that the localized-moment configuration has a lower energy in the DFT calculations presented here, which seems to contradict the conclusion that the distributed moment configuration is the relevant one for GaV 4 S 8 .However, the energy difference between the types of configurations is not large and it is possible that dynamical correlations may change the balance so that the distributed-moment case becomes lower in energy than the local moment configuration.In view of the calculated magnetic phase diagrams (Fig. 7) and their comparison with experiment, in particular, the observation of Néel skyrmions suggests strongly that the V 4 clusters are in the distributed-moment state.However, given the closeness of the different electronic configurations, and their distinctly different magnetic states, we speculate that the excitation spectra of GaV 4 S 8 should be particularly interesting, both from traditional electronand x-ray spectroscopic methods, as well as magnetic excitations, e.g., as provided by inelastic neutron scattering experiments.Experimental work is necessary to map out the complexities of the here proposed electronic and magnetic configurations of the V-based lacunar spinel GaV 4 S 8 .
FIG. 8. Electronic band structure of GaV4S8 calculated within density functional theory using LDA (a-c) and GGA (d-f) at different electronic correlation strengths characterized by the parameter U in the DFT+U scheme.For U = 0 eV the system is metallic but becomes semiconducting for U > 1 eV in LDA and for smaller U > 0 eV in GGA.For U = 2 eV, the results for the distributed-moment state are shown; this state is higher in energy than the localized-moment state, but both states show similar band gaps.Interestingly, the Dzyaloshinskii-Moriya interaction appears to change dramatically the temperaturedependence of the magnetization.This can be seen from the Monte-Carlo results obtained using only the Heisenberg interactions J ij (open symbols in Fig. 10a) compared to the simulations which include also the DM interaction D ij (filled symbols in Fig. 10a).The results suggest that the DMI enhances the long-range ferromagnetic order in a temperature range between 0 and 12 K, which is rather surprising considering that DMI usually tries to make the magnetic structure more non-collinear.In this respect, GaV 4 S 8 spinel contrasts many other known magnets and this aspect deserves further analysis which will be done in the future work.
Another important point is the difference between atomistic simulations based on inter-site (4-site model) or inter-cluster (1-site model) interactions, i.e. whether the individual sites or whole V 4 clusters are considered as the elementary magnetic units.From Fig. 10b, we can clearly see that neglecting the intracluster degrees of freedom by considering a 1-site spin model with effective interactions (eqn.( 2)) changes the behavior of the magnetization M (T ) significantly, as obtained in atomistic spin dynamics simulations, even when the DM interaction is not included.On the other hand, if the 1-site model includes the effective DM interaction, defined similarly to eqn.(2), we find that the long-range magnetic order is basically totally suppressed at any temperature above 1 K.In micromagnetic simulations, we do not see this drawback of the 1-site model.These findings suggest that the modelling of intracluster degrees of freedom in lacunar spinel GaV 4 S 8 is a subtle issue and has to be considered in greater detail in future studies.

Appendix C: Micromagnetic simulations
Based on the micromagnetic simulations for different temperatures and magnetic field strengths, we constructed the phase diagrams for GaV 4 S 8 with two different V 4 cluster states, shown in Fig. 7 in the main text.The original data points that we obtained in those mu-max3 simulations are depicted here in Fig. 11 using the same color code for convenience.It is worth noting that in our UppASD simulations using the multiscale module µASD 47 we observed skyrmion lattices in some conditions (example shown in Fig. 12).
For selected points in mumax3 phase diagrams, we have verified the effect of cell size along the z-direction, i.e. we compared the 2D and 3D simulation results.One example for external field of 250 mT is shown in Fig. 13 based on the simulation with 8 layers along the z-direction.It turns out that the phase diagrams do not change qualitatively, but the Curie temperature for the distributed-moment configuration increases up to 25 K, resulting in a shift of the boundary to the paramagnetic phase in Fig. 11a and agreeing better with the atomistic simulations in Fig. 4 in the main text.For the localizedmoment state (Fig. 11b), however, there is no visible shift of the Curie temperature in the 3D simulations.

FIG. 3 .
FIG. 3. Curie temperature estimated from classical atomistic Monte Carlo simulations based on first-principles Heisenberg Jij and Dzyaloshinskii-Moriya Dij intersite interactions calculated within GGA+U with variable correlation strength U .The dotted line represents the experimental Curie temperature.

FIG. 4 .
FIG. 4. Temperature-dependence of the magnetization M (T ) from Monte Carlo simulations based on magnetic interactions at U = 2 eV (solid lines), where both the localized and distributed moment states can be stabilized and compared.(30 × 30 × 30) supercell was used.The dotted line shows the M (T ) for the localized moment state at U = 1 eV.The shaded regions correspond to the magnetic susceptibility, the peak of which indicates the Curie temperature.

FIG. 5 .
FIG. 5. Schematic representation of the effective (clustercluster) DM vectors for nearest neighbors in GaV4S8 by arrows, where the arrow length is proportional to the DMI magnitude.The LSDA and LSDA+U results are shown in comparison.In both cases, the DM vectors follow the C3ν crystal symmetry.For the same bonds, the effective Heisenberg exchange (2) is given in Kelvin (positive value indicates a ferromagnetic interaction).

FIG. 6 .
FIG. 6. Micromagnetic simulation results (obtained with MuMax3 code) showing the annealed zero-temperature magnetic configurations: a) ferromagnetic domains with isolated skyrmions in external magnetic field B = 5 mT, b) zero-field spin-spiral state at T = 5 K, and c) Néel skyrmions with a size of 17 nm for B = 150 mT.The size of the simulation region is (512 × 512 × 1) nm in a,b) and (256 × 256 × 0.5) nm in c).The V4 clusters are in the distributed-moment state and the uniaxial anisotropy is set to K = 16 kJ/m 3 (other parameters are from the last row of TableII).The color code shows the out-of-plane direction of the magnetic moment vector in each point of space.

FIG. 7 .
FIG. 7. Theoretical phase diagrams of GaV4S8 for the a) distributed-and b) localized-moment V4 cluster configurations based on the micromagnetic simulations with MuMax3 code.Temperature and external magnetic field ranges are comparable to the experimentally studied ones.Different phases are shown: ferromagnetic (FM, orange), skyrmionic (SK, pink), cycloid+skyrmions (CY+SK, mint green), cycloid (CY, light blue) and paramagnetic (PM, yellow).The original data points for these phase diagrams are shown in Fig. 11 in the SI FIG. 9. Spin-resolved electronic band structure of GaV4S8 calculated within density functional theory in the GGA and GGA+U (U = 1 eV) formalisms.Red and blue curves correspond to spin-up and down d states projected on the V sites, and the V4 clusters are in the localized-moment configuration (Fig. 1b in the main text).

FIG. 10
FIG. 10. a) Temperature dependence of the magnetization M (T ) obtained from classical atomistic Monte Carlo simulations based on first-principles Heisenberg Jij and DM Dij intersite interactions (filled symbols) and Heisenberg interactions only (empty symbols), obtained within GGA+U = 1 eV.Different supercell sizes N 3 are used in the simulations.b) M (T ) for N = 30 obtained using the 4-site (solid lines) and 1-site (dashed lines) models, where in the latter case the intracluster degrees of freedom are not resolved.

FIG. 11 .
FIG.11.Original theoretical data points used for constructing the phase diagrams of GaV4S8 in Fig.7of the main text (the color code is the same).

FIG. 12 .
FIG. 12. Magnetic texture representing a) Néel-type skyrmion and b) a skyrmion lattice.The color bar represents the z component of the magnetization.The figure has been plotted by using the µASD software 47 .

TABLE II .
Calculated micromagnetic parameters (spin stiffness A in units of meV• ) are shown for comparison.