Ab initio exploration of short-pitch skyrmion materials: Role of orbital frustration

In recent years, the skyrmion lattice phase with a short lattice constant has attracted attention due to its high skyrmion density, making it a promising option for achieving high-density storage memory and for observing novel phenomena like the quantized topological Hall effect. Unlike conventional non-centrosymmetric systems where the Dzyaloshinsky-Moriya interaction plays a crucial role, the short pitch skyrmion phase requires a quadratic magnetic interaction $J(q)$ with a peak at finite-$Q$, and weak easy-axis magnetic anisotropy is also critical. Thus, conducting first-principles evaluations is essential for understanding the formation mechanism as well as for promoting the discovery of new skyrmion materials. In this {\it Perspective}, we focus on recent developments of the first-principles evaluations of these properties and apply them to the prototype systems Gd$T_2X_2$ and Eu$T_2X_2$, where $T$ denotes a transition metal and $X$ represents Si or Ge. In particular, based on the spin density functional theory with the Hubbard correction combined with the Liechtenstein method in the Wannier tight-binding model formalism, we first show that the Hubbard $U$ and Hund's coupling is essential to stabilize a skyrmion lattice state by enhancing the easy-axis anisotropy. We then discuss mechanisms of finite-$Q$ instability and show that competition among Gd-5$d$ orbitals determines whether ferromagnetism or a finite-$Q$ structure is favored in Gd$T_2$Si$_2$ with $T=$ Fe and Ru. Our systematic calculations reveal that GdRu$_2$$X_2$, GdOs$_2$$X_2$, and GdRe$_2X_2$ are promising, while GdAg$_2X_2$, GdAu$_2X_2$, and EuAg$_2X_2$ are possible candidates as the skyrmion host materials. Analysis based on a spin spiral calculation for the candidate materials is also presented.


I. INTRODUCTION
A magnetic skyrmion is a particle-like topological spin texture realized in condensed matter systems.Its magnetization distribution is characterized by vorticity, helicity, and a winding number, and the coupling to an itinerant electron provides a unique platform for exploring the physics under the emergent gauge field [1][2][3].The skyrmion is also considered a candidate element in future magnetoelectric devices due to its stability against an external perturbation and high controllability by an electric current [4][5][6][7][8].Particularly, a racetrack memory [9], a magnetic random access memory [10,11], and an artificial synapse for neuromorphic computing [12,13] are promising applications of the topologically stable magnetic object.
The magnetic skyrmions form a two-dimensional lattice in a ferromagnetic background.While an array of magnetic bubbles or skyrmions has been observed in a magnetic thin film since the 1970s [14,15], the first observation in a bulk magnet was accomplished in 2009 by the neutron scattering measurement for MnSi with B20 crystal structure [16].In the chiral magnet, an early theory suggested that a Dzyaloshinsky-Moriya (DM) interaction is the key to realizing the skyrmion lattice phase [17,18].Namely, the competition between ferromagnetic and DM interaction induces a non-collinear spin spiral structure, whose modulation pitch is determined by A/D with A be-ing magnetic stiffness and D being a DM constant.Once the spiral state is stabilized, a triple-Q state, which is identical to the triangular skyrmion lattice, is expected to appear under a magnetic field due to the momentum conservation of the quartic term in the Ginzburg-Landau functional.Today, there are many materials showing the DM-induced skyrmion phase in bulk, thin film, and heterostructures.In these systems, a typical modulation pitch is around 10-100nm, much larger than the lattice constant since D ≪ A is satisfied in most cases [19][20][21][22][23].
On the other hand, the magnetic skyrmion observed in the centrosymmetric crystals has recently attracted much attention, where the DM interaction should not play any role.The first observation was reported in Gd 2 PdSi 3 with Gd forming a two-dimensional triangular lattice network [24,25], followed by Gd 3 Ru 4 Al 12 with a Kagome lattice [26], GdRu 2 Si 2 [27], EuAl 4 [28][29][30][31], and EuGa 4 [30,32] with a square lattice network.A remarkable feature of the skyrmion lattice in centrosymmetric systems is its relatively shorter modulation pitch, typically given by a few nanometers.Since this does not originate from the DM interaction, its formation mechanism has been intensively studied recently.In terms of the shortness of the modulation pitch, the skyrmion lattices observed in EuPtSi [33,34] and Mn 1−x Fe x Ge [35,36] might be classified into the same category as the centrosymmetric systems even though their crystal structures do not have the space inversion symmetry.Since a short-pitch skyrmion lattice is favorable for realizing high-density memory bits in a storage device, revealing the key quantity to determine the modulation pitch will promote further materials design and engineering of the arXiv:2404.02428v1[cond-mat.str-el]3 Apr 2024 skyrmion phase.Furthermore, the high skyrmion density of the short-pitch skyrmion is of great theoretical importance, as it provides a unique platform to observe fascinating phenomena, such as a quantized topological Hall effect [37] and a crossover behavior of the topological Hall effect where both the real-space and momentumspace Berry curvature play a crucial role [38].
To date, there have been a lot of theoretical works on the formation of a short-pitch skyrmion lattice.They include calculations based on a classical spin model, an itinerant electron model, and an effective Ginzburg-Landau functional.These theories enable us to grasp a physical insight into the skyrmion lattice formation and, indeed, the importance of the magnetic frustrations [39,40], compass anisotropy of the magnetic interactions [41,42], higher-order spin interactions [43,44], and Fermi surface nesting [45,46] have been recognized as an essential property to stabilize the skyrmion lattice phase (for more details, see Ref. [47]).However, these model calculations are usually performed with a given modulation vector and include adjustable parameters.Thus, they have no predictive power for the modulation pitch itself in real materials.On the other hand, in spite of the expensive numerical cost, one can investigate the stability of finite-Q spin spiral structures from first-principles, which is helpful to see the origin of the modulation pitch and its material dependence.Indeed, the calculations for GdRu 2 Si 2 [48,49], Gd 2 PdSi 3 [48,49], and MnGe [50,51] have successfully unveiled microscopic mechanisms in stabilizing their short-pitch skyrmion phase.Thus, broadening the target compounds along this line is of great importance, helping us understand the material dependence more deeply and promoting the discovery of new skyrmion materials.
Based on the theory for short-pitch skrymion phase, a promising system for observing skyrmion phases should possess weak easy-axis anisotropy and spin-spin interactions favoring finite-Q modulations [47].Since these effects on skyrmion formation are well-established, in this Perspective, we aim to highlight recent progress in firstprinciples evaluations of these properties and offer predictive insights to guide experimental efforts in discovering new skyrmion materials.To achieve this, based on spin density functional theory (SDFT) and SDFT with the Hubbard correction U (SDFT+U ), we perform systematic first-principles calculations for GdT 2 X 2 and EuT 2 X 2 with T being a transition metal element and X being Si or Ge.Since the largest number of skyrmion compounds have been found in this ThCr 2 Si 2 -type crystal structure, these compounds would be ideal target materials.First, we show that the inclusion of U and Hund's coupling J is essential to obtain easy-axis anisotropy, which is necessary for stabilizing a skyrmion lattice than a spin spiral state.Then, we evaluate magnetic interactions J(q) based on the so-called Liechtenstein method with the Wannier tight-binding model formalism.The calculated J(q) is in good agreement with E(q) directly evaluated by the spin spiral calculations within SDFT+U , support-ing the validity of the present Liechtenstein calculations.Based on the results, we argue that the finite-Q structure is determined not only by the shape of the Fermi surface but also by the details of the electronic structure.Especially, we show that competition among the Gd-5d orbital contributions determines whether a ferromagnetism or finite-Q structure is favored in GdT 2 Si 2 with T = Fe and Ru.According to our systematic calculations, and EuAg 2 X 2 with X = Si and Ge are possible candidates for showing the skyrmion phase.The present study provides a firm ground to discover and design the short-pitch skyrmion phase from systematic calculations, and the extension to the other crystal structures will be a promising future development.

II. METHODS
In this section, we summarize the methods used in this Perspective.First, we explain the evaluation of J(q) based on the Liechtenstein method.Although this method was first formulated in the multiple scattering theory with the Green's function techniques, it has been applied to the SDFT Hamiltonian with various spatially localized bases, including the Wannier functions [52,53].Here, we show the formulation following Ref.[53].We begin with the following tight-binding Hamiltonian, where indices 1 and 2 run over all degrees of freedom that specify the Wannier functions, including atomic sites i, orbitals ℓ, and spins σ. c † 1 (c 1 ) represents an electron creation (annihilation) operator in this basis.t 12 and v 12 denote spin-independent and spin-dependent hopping integral matrices, respectively.These parameters are extracted from the SDFT/SDFT+U calculation with ferromagnetic reference states through the Wannier construction process [54,55].
Based on the Hamiltonian (1), we can evaluate the magnetic interaction J ij in the classical Heisenberg model as follows, Here, Tr nℓσ = T n Tr ℓσ , where Tr ℓσ denotes the trace over the orbital and spin spaces.ω n = πT (2n+1) denotes the Matsubara frequency, and T is the temperature.The Green's function G(ω n ) and the magnetic perturbation matrix at i site Σ i are defined by, where we approximate Σ i 12 as a local quantity.σ µ (µ = x, y, z) is the Pauli matrix, and ṽi1ℓ1,i2ℓ2 is de- is defined as a sub-matrix of G 12 (ω n ) having (i, j) sites component.The evaluation of the Matsubara summation in Eq. ( 2) is performed by the sparse sampling technique with the intermediate representation [56,57].More details about the implementation are found in Ref. [53].
Based on J ij evaluated by Eq. ( 2), its Fourier transform J(q) is obtained as follows: where R i denotes the real space coordinate of the site i including the relative coordinate of the sublattice.Here, we have introduced J 0 , the on-site magnetic interaction, to guarantee the absence of the self-interaction term in the Heisenberg model.Note that, although J ij evaluated by Eq. ( 2) will be finite for non-magnetic sites such as a transition metal T and X = Si, Ge in GdT 2 X 2 and EuT 2 X 2 , the i, j summation in Eq. ( 5) should be taken only for the magnetic site, namely, Gd or Eu site in our cases.Indeed, J ij connecting non-magnetic elements can take a large value when both the density of states (DOS) and exchange splitting near the Fermi level are finite, although the mapping into the classical spin model for these orbitals is totally inadequate.This procedure can remove this artificial contribution in the evaluation of J(q), which is essential for the SDFT+U calculations.
Based on the obtained J(q), we can see whether the system favors the finite-Q modulation.
On the other hand, we can evaluate the energy of the spin spiral states with the modulation vector Q, denoted by E(q), within SDFT and SDFT+U .The calculation is based on the generalized Bloch theorem, in which the energy of the spin spiral state is evaluated with specific twisted boundary conditions for the Bloch spinors [58].A clear advantage of this method is that the calculated E(q) is exact within the SDFT/SDFT+U level, and no assumption of the mapping to the classical spin model is introduced.Moreover, it is not necessary to expand the unit cell even when calculating the finite-Q structure, which makes the calculation with small Q tractable.However, despite these advantages, it requires a calculation for every Q, and thus, is much more expensive than the Liechtenstein method if one wants to see the detailed structures of E(q).Thus, we use these two methods in a complementary manner.

III. CALCULATION DETAILS
For the evaluations of J(q) and E(q), we use Vienna ab initio simulation package (VASP) [59] for the electronic structure calculations.Here, the exchangecorrelation functional proposed by Perdew, Burke, and Ernzerhof [60], and pseudopotentials with the projector augmented wave (PAW) basis [61,62] are used.The spin orbit interaction is neglected except for the magnetocrystalline anisotropy calculations.The structure optimization is performed for all GdT 2 X 2 and EuT 2 X 2 compounds by the SDFT+U method assuming ferromagnetic states.The convergence criteria of the structure optimization is set to 10 −5 eV and the corresponding electronic self-consistent loop is 10 −6 eV with the accurate precision mode.Here, we use a 16 × 16 × 16 Monkhorst-Pack k-mesh, and a cutoff of the plane wave basis set is set to be twice the recommended maximum value of the cutoff.For the magnetocrystalline anisotropy calculations, we employ 10 −8 eV as the convergence criteria for electronic calculations and 21 × 21 × 21 Monkhorst-Pack k-mesh.The cutoff is set to be twice the recommended value.For the energy calculations of the spin spiral states, we employ 10 −8 eV as the convergence criteria for electronic calculations, 16 × 16 × 16 Monkhorst-Pack k-mesh, and the cutoff is 2.5 times as large as the recommended values.
For constructing the Wannier tight-binding model, we use 156 Bloch states evaluated on 10 × 10 × 10 uniform kgrid in the disentanglement procedure.The trial orbitals are set to the Gd/Eu-d and f orbitals, transition metal T -s, p, and d orbitals, and X = Si/Ge-s and p orbitals.The outer window is set as it includes all 156 Bloch states, and the inner window is set as it covers from the bottom of the bands to the bands at 4 eV above the Fermi level.The obtained tight-binding model consists of 38 orbitals per spin.Based on the obtained tight-binding model, J ij and J(q) are evaluated by Eqs. ( 2) and ( 5), respectively.Here, we set the temperature T = 58 K and employ 48 × 48 × 48 uniform k-grid.

A. Magnetocrystalline anisotropy
According to the theory of skyrmion lattice formation [47], the system should possess two properties to realize a skyrmion lattice phase in a broad region of parameter space.Namely, it is necessary that the spin in-teraction is compatible with a finite-Q modulation, and the magnetocrystalline anisotropy is the easy axis type as well as sufficiently weak.It should be noted that, with the exception of Gd 3+ and Eu 2+ , rare earth elements possess a finite orbital moment with strong spin-orbit coupling, leading to strong magnetocrystalline anisotropy in most cases.As a result, compounds containing these elements are unlikely to exhibit the weak magnetocrystalline anisotropy required for the formation of skyrmion lattice phases.Therefore, in this Perspective, we only consider Gd 3+ -and Eu 2+ -based compounds, where the orbital angular momentum of their 4f -orbitals are almost quenched in the ionic ground states.On the other hand, there is no general trend for these elements, whether the anisotropy is the easy axis or easy plane type.Thus, it is crucial to see whether the system shows the correct magnetic anisotropy before discussing the modulation vectors based on the electronic structure.
As an example, let us first see in Fig. 1 the calculated magnetocrystalline anisotropy, for the skyrmion compound GdRu 2 Si 2 and a reference compound GdRu 2 Ge 2 [63].We can see that the energy scale of the magnetic anisotropy is quite small, |∆E| < 5 meV per formula unit, as expected.However, the sign of ∆E is negative in the SDFT results, indicating that a helical magnetic structure is expected to have a lower energy than the skyrmion state.This situation drastically changes by including the Hubbard correction U .In the case of Gd 3+ , since U enhances the energy difference between the lowest energy state with L = 0 and the first excited state with L ̸ = 0, we can expect that the magnetocrystalline anisotropy decreases by increasing the value of U , which is consistent with the results in Fig. 1.Notably, ∆E becomes small but positive with U = 6.7 eV even without Hund's coupling correction J.
In addition, we can also see that the anisotropy becomes larger positive with increasing the strength of J, which is compatible with the observed skyrmion lattice phases.This result clearly indicates that the inclusion of U and J is essential to obtain the correct magnetocrystalline anisotropy in the Gd and Eu-based skyrmion materials.Figures 2(a-d) summarize the results of ∆E in all Gd and Eu-based 122 compounds that we discuss for this Perspective.We can see that, for most transition metal T except for Ta, the SDFT+U results show the same or lower |∆E| than SDFT, which is similar to the cases of GdRu 2 Si 2 and GdRu 2 Ge 2 .Although the values of ∆E highly depend on the host f -element and Eu-based compounds tend to have considerable negative ∆E in SDFT, these trends are strongly suppressed in SDFT+U .In SDFT+U , ∆E is no longer sensitive to the transition metal T and the group 14 element X = Si and Ge.Also, the amplitude |∆E| is sufficiently small and |∆E| < 3 meV/f.u. is satisfied in all compounds.It is worth noting that only one (seven) in GdT 2 X 2 (EuT 2 X 2 ) among all 48 compounds considered show negative ∆E, implying that the weak easy-axis anisotropy is an intrinsic property inherent in this ThCr 2 Si 2 -type crystal structure.

B. Fermi surface and magnetic instability
Next, let us move on to the modulation vector of GdT 2 X 2 and EuT 2 X 2 .According to the previous discussion, it is essential to include the local Coulomb repulsion U and Hund's coupling J to reproduce the correct magnetocrystalline anisotropy observed in the experiments.Thus, in the present study, we discuss the modulation vector Q based on the electronic structures in the SDFT+U calculations.Apparently, the effect of U and J on Q should not be significant if Q is determined only by the shape of the Fermi surface, which is often supposed in the RKKY scenario for the skyrmion lattice formation.Indeed, since the magnetic moments of Gd 3+ and Eu 2+ are fully polarized even in the absence of U and J, the Fermi surface is expected to be insensitive to the values of U and J.We can directly confirm that magnetic moments are fully polarized in the ferromagnetic states in most compounds except for EuMn 2 Si 2 , EuZr 2 Si 2 , EuZr 2 Ge 2 , EuTc 2 Si 2 , EuTc 2 Ge 2 , EuHf 2 Ge 2 , EuTa 2 Si 2 and EuRe 2 Ge 2 .Thus, although the electronic structure above the Fermi energy should be modified to some extent because of the change of the exchange splitting, this will not affect the modulation vector Q in the skyrmion phase.
The insensitivity of the Fermi surface against U and J can be seen directly by the band structure calculations.For example, we show the band structures near the Fermi levels of GdRu 2 Si 2 and GdFe 2 Si 2 in Figs.3(a-d).Here, the majority and minority bands calculated with SDFT and SDFT+U are plotted, respectively.We can see from the figures that, although they differ in detail especially above the Fermi levels, U and J do not change the Fermi surface, which means that SDFT and SDFT+U results share the same Fermi surface instability.Note that this trend is valid for all compounds when SDFT correctly captures the ionic state of Gd 3+ and Eu 2+ .
In the continuum model, it is well-known that the Lindhard function describes the Fermi surface instability.In a real material with many bands, the nesting function, plays the same role.Here, n and k denote the band index and crystal momentum, respectively.ε nk is the eigenvalues of the Bloch states having n and k, and f (ε) is the Fermi distribution function.The nesting function is sometimes considered to describe not only the Fermi surface instability but also the spin instability, such as spin-density-wave and helical Q magnetic structures.However, as is known in the fluctuation theory for multi-orbital systems, the orbital character of the bands is essential to describe the fluctuation.In this case, it is stand for the calculations with U = 6.7 eV and J = 0.7 eV.Note that the data for EuMn2Si2, EuZr2Si2, EuZr2Ge2, EuTc2Si2, EuTc2Ge2, EuHf2Ge2, EuTa2Si2 and EuRe2Ge2 in SDFT, and GdCr2Ge2 in SDFT+U are not shown.This is because the fully polarized ferromagnetic solutions along the x and/or z axis are unstable for these materials.
FIG. 3. Band structures of (a) majority spin of GdRu2Si2, (b) minority spin of GdRu2Si2, (c) majority spin of GdFe2Si2 and (d) minority spin of GdFe2Si2.Solid lines correspond to SDFT and dashed lines correspond to SDFT+U results with U = 6.7 eV and J = 0.7 eV.
necessary to consider the (zeroth-order) spin correlation function χ µν (q) defined as follows, e iϕ q 13 (R) s µ 13 s ν 24 χ 13,24 (R), ( 8) Here, the phase ϕ q 12 (R) is defined by ϕ q 12 (R) = q • (R i3 − R i1 ), and ω n and ω q denote the Fermionic and Bosonic Matsubara frequency, respectively.Since Eq. ( 9) leads to the nesting function with a modification due to the transformation from local orbital to band basis, the qdependence of χ µν (q) is mainly determined by Eq. ( 7).However, it should be noted that the basis transformation gives an additional source of the q-dependence, and χ µν (q) is sensitive not only to the shape of the Fermi surface but also to its orbital character.Moreover, in the strongly correlated electron systems, Eqs. ( 8) and ( 9) are not sufficient to describe the correct spin fluctuation.In this case, the effect of U must be taken more accurately using a diagrammatic technique, which is practically impossible in most cases.In that sense, the Liechtenstein method presented provides a simple but reliable approximation in taking the correlation effects.Indeed, it is known that this method successfully reproduces the result of t/U expansion in the strong correlation limit and reproduces Eqs. ( 8) and (9) in the weak correlation limit.It is worth noting that in the intermediate and strong U regime, not only the Fermi surface but also the electronic structure far from the Fermi level can affect the spin instability.
Figures 4(a) and (b) show J(q) for GdRu 2 Si 2 and GdFe 2 Si 2 calculated based on the Liechtenstein formula, Eq. ( 5).The result in Fig. 4(a) shows good agreement with that in Ref. [49] but is slightly different from that in Ref. [48] due to an erroneous treatment of the phase factor in Eq. (8).According to the results, since the peak position Q in J(q) is the signature of magnetic instability with the modulation vector Q, we can say that GdRu 2 Si 2 favors the spin spiral state with Q ∼ (0.25, 0, 0) both in SDFT and SDFT+U .On the other hand, in GdFe 2 Si 2 , finite-Q instability is seen only in SDFT while SDFT and SDFT+U give the same Fermi surface (see Figs. 3(c) and (d)).This result indicates that the finite-Q structure is determined not only by the shape of the Fermi surface but also by other factors, such as an orbital character and electronic structure far from the Fermi level.The ferromagnetic instability observed in the SDFT+U calculation for GdFe 2 Si 2 is also found in the spin spiral calculation (which will be discussed in Section IV F).However, it seems to be inconsistent with the observed spin flop transition in GdFe 2 Si 2 [64], which implies that we may need to treat the correlation effect more precisely to understand the magnetisn in GdFe 2 Si 2 .On the other hand, GdRu 2 Si 2 has a peak at Q ∼ (0.25, 0, 0) both in SDFT+U and SDFT calculation, similar to the previous results based on SDFT (without U ) [48,49], which is consistent with the experiment.
Note that, in principle, it is necessary to perform a spin model simulation to find out a correct magnetic ground state, especially for multiple-Q states like a skyrmion lattice phase.For a spin model of GdRu 2 Si 2 , an analysis based on the Landau-Lifshitz-Gilbert (LLG) equation was performed in Ref. [49], using the following model for the internal energy, where K is a magnetocrystalline anisotropy constant, and B is a external magnetic field applied along the c axis.They solved the LLG equations for Eq. ( 10) with J ij evaluated from first-principles, which is consistent with Fig. 4(a), and obtained a magnetic phase diagram in terms of K and B, as shown in Fig. 5.According to Fig. 5, the skyrmion lattice phase appears when the easyaxis anisotropy K is in the range from 0.2 to 0.4 meV.
On the other hand, as discussed in the previous sections, the SDFT+U calculation reproduces the easy-axis anisotropy of GdRu 2 Si 2 .The value of K in GdRu 2 Si 2 is ∼ 0.43 meV, which is in good agreement with the spin model simulation Note that the significance of weak easyaxis anisotropy can be straightforwardly understood: if the anisotropy is easy-plane type, a helical single-Q state is favored over multiple-Q states like a skyrmion phase, and if the easy-axis anisotropy is too strong, Ising-type magnetic order is favored over finite-Q structures.Thus, we conclude that the SDFT+U result in GdRu 2 Si 2 is fully consistent with the observation of the skyrmion lattice phase, and the present method correctly captures the essential physics behind the magnetic structure.

D. Orbital decomposition of J(q) calculated based on the Liechtenstein method
One of the advantages of evaluating J(q) based on the tight-binding formalism is that it is easy to discuss the microscopic origin of the magnetic interaction [48,65].Although J(q) was decomposed into the Gd-4f and 5d components and a possible frustration mechanism between these two was discussed in Ref. [48], one may think that the situation will be different in a real material since their calculations were based on the SDFT electronic structures.Indeed, with the SDFT+U results, we can easily show that the contribution from Gd-4f orbital is strongly suppressed and never affect the modulation vector of the helical states, as shown in Fig. 6(a).Naively, this indicates the absence of conventional RKKY interac- tions in this compound since it is included in the Gd-4f contribution as discussed in Ref. [48].On the other hand, as seen in the previous subsections, the orbital character and the electronic structures not at the Fermi level are essential factors in determining the finite-Q instability.
Since the q dependence of J(q) mainly comes from the Gd-5d manifold according to Fig. 6(a), it would be interesting to decompose J(q) into each 5d orbital based on the SDFT+U electronic structures.Figure 6(b) and (c) show the diagonal and off-diagonal contributions to J(q) in the Gd-5d manifold, respectively.Among the five 5d orbitals, the contributions from the d z 2 and d xy orbitals are negligibly small, and the q dependence of J(q) originates from the remaining three orbitals.Among the three, we can see that only the diagonal contribution from the d x 2 −y 2 orbital has a peak structure at the Γ point, favoring the ferromagnetic ground state.In contrast, the other diagonal and off-diagonal ones show the peak at Q ∼ (0.25, 0, 0) or almost flat along the Γ-X line.These results indicate that the orbital frustration exists in the 5d orbital manifold, and (0, 0, 1  2 ) (0, 0, 0) ( 1 2 , 0, 0) T = Ru total 44 others T = Fe total 44 others FIG. 7. Orbital decomposed J(q) calculated by the Liechtenstein method.The total (dotted line with open square), diagonal d x 2 −x 2 orbital (solid line with filled square), and the other contribution (dotted like with filled circle) are shown.The pink and cyan lines correspond to the results for GdRu2Si2 and GdFe2Si2, respectively.
whether the system favors the ferromagnetic or finite-Q spin spiral state depends on the competition between the contributions from d x 2 −y 2 and the other orbitals.Notably, this situation is realized not only in GdRu 2 Si 2 but also in GdFe 2 Si 2 as is shown in Fig. 7.The figure shows why these two compounds favor different magnetic structures though they share almost the same Fermi surface.Namely, the contribution from the d x 2 −y 2 orbital relative to that from the other d orbitals is larger in GdFe 2 Si 2 than GdRu 2 Si 2 .This is mainly due to the suppression of the finite-Q peak in the contribution from the other d orbitals in GdFe 2 Si 2 .
Here, let us look into the detail of the microscopic origin of J(q) based on the following simple model Hamiltonian H = H d + H f + H df , where, Here, H d , H f and H df represent terms for the Gd-5d orbitals, 4f orbitals and their hybridization terms, respectively.Our calculations for GdT 2 X 2 show that, in the presence of U , the 5d contribution always overcomes that of 4f since the latter is strongly suppressed.This clearly indicates that both the super-exchange interaction J f ex ∼ t 2 f /U f and the conventional RKKY interaction J RKKY ∼ (V 2 /U f ) 2 χ d , where χ d is the spin susceptibility in the 5d orbitals, are negligibly small in these compounds.Conversely, the strong magnetic interaction from the 5d orbitals can be easily understood since they have large DOS near the Fermi level, and have a finite exchange splitting due to the magnetic ordering although it is not as large as that of the 4f orbitals.Here, we note that there are two possibilities for the origin of the exchange splitting in the 5d states.One is the Coulomb interaction inherent in the 5d orbitals, i.e., the second term in eq. ( 11), in which case the spin instability is described purely by the 5d orbitals.The other is the magnetic coupling between the 5d and 4f orbitals, i.e., the second term in eq. ( 13), which one may call "generalized" RKKY interaction since it comes from the coupling between conduction bands and local spins [66].It should be noted that in SDFT+U , the correlation effect is included in the exchange-correlation functional, and the many-body problem turns into an effective onebody calculation.Then, we can not distinguish these two contributions from the resulting exchange splittings.Investigating the microscopic origin of J(q) more precisely requires more elaborated calculations, such as a calculation combined with dynamical mean-field approximation [67].This direction of development is an interesting and important future task in this field.

E. Systematic evaluation of J(q) based on the Liechtenstein method
Let us move on to the results of J(q) calculated with SDFT+U for GdT 2 X 2 and EuT 2 X 2 in Fig. 8 and Fig. 9, respectively.From these figures, we see that the general trend of J(q) is not sensitive to the choice of group 14 elements X = Si or Ge.On the other hand, the choices of different f -elements and different periods (i.e., 3d, 4d, or 5d) of the transition metals T dramatically change the structure of J(q).For example, for the GdT 2 X 2 compounds, the T dependence of group 5 (T = V, Nb, Ta), 8 (T = Fe, Ru, Os), and 9 (T = Co, Rh, Ir) are not so significant.On the other hand, for EuT 2 X 2 , the T dependence is more noticeable.These dependencies may come from the chemical pressure effects on the lattice structures or the spreads and local energy levels of the d-orbitals.
According to Figs. 8 and 9, the most promising candidates showing skyrmion lattice phase would be compounds with T being the group 8 transition metal.In particular, GdRu 2 Si 2 , GdRu 2 Ge 2 , GdOs 2 Si 2 , GdOs 2 Ge 2 , and EuOs 2 Si 2 have the highest peaks at finite Q along the Γ-X line.Together with the small easy-axis anisotropy shown in Figs.2(a-d), they should show the skyrmion lattice phase, and indeed, it is already found in GdRu 2 Si 2 .Except for these compounds, we can also see that GdW 2 X 2 , GdRe 2 X 2 , GdAg 2 X 2 , GdAu 2 X 2 , EuCo 2 X 2 , and EuAg 2 X 2 with X = Si and Ge have the highest peaks along the Γ-X line.However, the peaks in GdW 2 X 2 and EuCo 2 X 2 may be too small to stabilize the magnetic order.

F. Modulation vector based on spin spiral calculations
Up to now, we have discussed the magnetic interaction J(q) and the corresponding modulation vectors Q based on the Liechtenstein method.The calculations are performed for the Wannier tight-binding models derived from first-principles for the ferromagnetic ground states.Thus, in principle, the reliability of the q-dependence far from the Γ point, as well as the validity of the mapping to the classical spin model cannot be justified.Here, we present spin spiral calculations for the materials expected to be a good candidate for showing the skyrmion phase.Since the calculation is exact within the SDFT+U level, the results can be used in a way complementary to that of the Liechtenstein method.
Figure 10(a) shows the q-dependence of the total energies for GdT 2 X 2 with T = Fe, Ru, Os.The calculations are performed by the spin spiral calculations based on SDFT+U with U = 6.7 and J = 0.7 eV.Since the minimum position of E(q) represents the most stable modulation vector, we can see that GdRu 2 Si 2 , GdRu 2 Ge 2 , and GdOs 2 Si 2 favor the finite-Q state with Q ∼ (0.20, 0, 0)-(0.25,0, 0).On the other hand, GdFe 2 Si 2 and GdFe 2 Ge 2 have the minimum at the Γ point, indicating that the ferromagnetic state is the most stable.These features are in good agreement with the calculations based on the Liechtenstein method.Although the GdOs 2 Ge 2 does not show the finite-Q instability in contrast with the Liechtenstein result, the peak of J(q) in GdOs 2 Si 2 is more fragile than the others as is shown in Fig. 8(e).Thus, the disagreement with the Liechtenstein method is not serious.For the comparison, we also show E(q) of GdT 2 X 2 with T = Fe, Ru, Os and X = Si, Ge based on SDFT in Fig. 10(b).We can see that the finite-Q instability becomes stronger in all cases than in SDFT+U .Here, GdFe 2 Si 2 shows the finite-Q instability in SDFT while it does not in SDFT+U , which again agrees well with the Liechtenstein calculations.In Fig. 10(c), we show E(q) for the other candidates for showing the skyrmion phase predicted by the Liechtenstein calculations, namely, GdT 2 X 2 with T = Re, Ag, Au and X = Si, Ge, EuAg 2 Si 2 and EuAg 2 Ge 2 .We can see that GdRe 2 Si 2 , and GdRe 2 Ge 2 show nearly flat q dependence near the Γ point, and thus, it is possible that some subtle perturbations select a finite-Q state.On the other hand, GdAg 2 X 2 , GdAu 2 X 2 , and EuAg 2 X 2 have a broad peak at Q ∼ ( 1 2 , 0, 0) corresponding to the 90degree rotating structure with the nearest neighboring spins, which seems to be too large to realize a skyrmion phase.However, due to its broad nature, we may still have a chance to achieve a skyrmion phase with smaller Q in experiments by using, for example, chemical substitution and external pressure.We leave the possible magnetic structures favored in these compounds as a future study.8. Spin interactions J(q) in GdT2X2 calculated by the Liechtenstein method based on the SDFT+U electronic structures.The pink, green, and cyan lines with square, circle, and triangle points correspond to 3d, 4d, and 5d elements as T , respectively.The solid (dashed) lines with filled (open) symbols stand for X = Si (Ge).Note that the high symmetry points Z, Γ, and X correspond to (0, 0, 2π c ), (0, 0, 0) and ( 2π a , 0, 0) in the cartesian frame of the reciprocal lattice space, respectively.

V. CONCLUSION AND OUTLOOK
In this Perspective, we present systematic firstprinciples calculations for GdT 2 X 2 and EuT 2 X 2 with T being a transition metal element and X being Si or Ge.From the magnetocrystalline anisotropy calculations based on SDFT and SDFT+U , we show that the inclusion of Coulomb interaction U and Hund's coupling J   is essential to obtain the easy-axis anisotropy.Then, based on the SDFT+U electronic structures, we evaluate magnetic interactions J(q) based on the Liechtenstein method and show that the obtained J(q) is in good agreement with E(q) by the spin spiral calcula-tions.Our calculations indicate that the finite-Q structure is determined not only by the Fermi surface topology but also by the details of the electronic structure, and the competition of each Gd-5d orbital contribution determines whether a ferromagnetic spin configuration or finite-Q structure is favored in GdT 2 Si 2 with T = Fe and Ru.According to our calculations, GdRu 2 X 2 , GdOs 2 X 2 , and GdRe 2 X 2 are promising candidates, while GdAg 2 X 2 , GdAu 2 X 2 , and EuAg 2 X 2 are possible candidates for showing the skyrmion lattice phase.Since the systematic calculations based on the Liechtenstein method is shown to be helpful in evaluating the finite-Q structure, the extension to other crystal structure would be a possible direction to discover and engineer new skyrmion compounds.On the other hand, developing new methods to evaluate transport properties in the short pitch skyrmion phase is another crucial direction for its practical applications.Since most of the novel phenomena associated with its high skyrmion density have been studied based on simple models, further development is necessary to achieve a better understanding of material dependence and quantitative evaluations.

2 FIG. 2 .
FIG.2.Magnetocrystalline anisotropy of (a) GdT2Si2, (b) GdT2Ge2, (c) EuT2Si2 and (d) EuT2Ge2.The SDFT+U results stand for the calculations with U = 6.7 eV and J = 0.7 eV.Note that the data for EuMn2Si2, EuZr2Si2, EuZr2Ge2, EuTc2Si2, EuTc2Ge2, EuHf2Ge2, EuTa2Si2 and EuRe2Ge2 in SDFT, and GdCr2Ge2 in SDFT+U are not shown.This is because the fully polarized ferromagnetic solutions along the x and/or z axis are unstable for these materials.

FIG. 4 .
FIG. 4. Modulation vector Q dependence of the spin exchange interaction J(q) in (a) GdRu2Si2 and (b) GdFe2Si2 calculated based on the Liechtenstein formula.A pink (cyan) line with open (filled) circle corresponds to the SDFT (SDFT+U with U = 6.7 eV and J = 0.7 eV) result.

FIG. 6 .
FIG.6.Orbital decomposed J(q) calculated by the Liechtenstein method.(a) The contributions from the Gd-5d and 4f orbitals.(b) The diagonal contribution to in the Gd-5d manifold.Here, the indices 1 to 5 correspond to the d z 2 , dxz, dyz, d x 2 −y 2 , and dxy orbitals, respectively.(c) The off-diagonal contributions in the Gd-5d manifold.
FIG.8.Spin interactions J(q) in GdT2X2 calculated by the Liechtenstein method based on the SDFT+U electronic structures.The pink, green, and cyan lines with square, circle, and triangle points correspond to 3d, 4d, and 5d elements as T , respectively.The solid (dashed) lines with filled (open) symbols stand for X = Si (Ge).Note that the high symmetry points Z, Γ, and X correspond to (0, 0, 2π c ), (0, 0, 0) and ( 2π a , 0, 0) in the cartesian frame of the reciprocal lattice space, respectively.

FIG. 9 .
FIG. 9. Spin interactions J(q) in EuT2X2 calculated by the Liechtenstein method based on the SDFT+U electronic structures.Color, line and point types are the same as Fig.8.

FIG. 10 .
FIG.10.Modulation vector q dependence of the ground state energy E(q) calculated based on the Frozen magnon method.(a) The results based on the SDFT+U with U = 6.7 eV and J = 0.7 eV for GdT2X2 with T = Fe, Ru, Os and X = Si, Ge.(b) The results based on the SDFT for GdT2X2 with T = Fe, Ru, Os and X = Si, Ge.(c) The results based on the SDFT+U with U = 6.7 eV and J = 0.7 eV for GdT2X2 with T = Re, Ag, Au and X = Si, Ge.