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

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 interaction 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³⁺ and Eu²⁺, 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³⁺- and Eu²⁺-based compounds, where the orbital angular momentum of their 4$f$-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₂Si₂ and a reference compound GdRu₂Ge₂ [63]. We can see that the energy scale of the magnetic anisotropy is quite small, $|\Delta E|<5$ meV per formula unit, as expected. However, the sign of $\Delta 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³⁺, since $U$ enhances the energy difference between the lowest energy state with $L=0$ and the first excited state with $L\neq 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, $\Delta 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 $\Delta 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 $|\Delta E|$ than SDFT, which is similar to the cases of GdRu₂Si₂ and GdRu₂Ge₂. Although the values of $\Delta E$ highly depend on the host $f$-element and Eu-based compounds tend to have considerable negative $\Delta E$ in SDFT, these trends are strongly suppressed in SDFT+$U$. In SDFT+$U$, $\Delta E$ is no longer sensitive to the transition metal $T$ and the group 14 element $X=$ Si and Ge. Also, the amplitude $|\Delta E|$ is sufficiently small and $|\Delta E|<3$ meV/f.u. is satisfied in all compounds. It is worth noting that only one (seven) in Gd$T_{2}X_{2}$ (Eu$T_{2}X_{2}$) among all 48 compounds considered show negative $\Delta E$, implying that the weak easy-axis anisotropy is an intrinsic property inherent in this ThCr₂Si₂-type crystal structure.

Next, let us move on to the modulation vector of Gd$T_{2}X_{2}$ and Eu$T_{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 $\bm{Q}$ based on the electronic structures in the SDFT+$U$ calculations. Apparently, the effect of $U$ and $J$ on $\bm{Q}$ should not be significant if $\bm{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³⁺ and Eu²⁺ 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₂Si₂, EuZr₂Si₂, EuZr₂Ge₂, EuTc₂Si₂, EuTc₂Ge₂, EuHf₂Ge₂, EuTa₂Si₂ and EuRe₂Ge₂. 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 $\bm{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₂Si₂ and GdFe₂Si₂ 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³⁺ and Eu²⁺.

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 $\bm{k}$ denote the band index and crystal momentum, respectively. $\varepsilon_{n{\bm{k}}}$ is the eigenvalues of the Bloch states having $n$ and $\bm{k}$, and $f(\varepsilon)$ 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 $\bm{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 necessary to consider the (zeroth-order) spin correlation function $\chi_{\mu\nu}(\bm{q})$ defined as follows,

Here, the phase $\phi_{12}^{\bm{q}}(\bm{R})$ is defined by $\phi_{12}^{\bm{q}}(\bm{R})=\bm{q}\cdot(\bm{R}_{i_{3}}-\bm{R}_{i_{1}})$, and $\omega_{n}$ and $\omega_{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 $\bm{q}$-dependence of $\chi_{\mu\nu}({\bm{q}})$ is mainly determined by Eq. (7). However, it should be noted that the basis transformation gives an additional source of the $\bm{q}$-dependence, and $\chi_{\mu\nu}({\bm{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(\bm{q})$ for GdRu₂Si₂ and GdFe₂Si₂ 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 $\bm{Q}$ in $J(\bm{q})$ is the signature of magnetic instability with the modulation vector $\bm{Q}$, we can say that GdRu₂Si₂ favors the spin spiral state with $\bm{Q}\sim(0.25,0,0)$ both in SDFT and SDFT+$U$. On the other hand, in GdFe₂Si₂, finite-$\bm{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-$\bm{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₂Si₂ is also found in the spin spiral calculation (which will be discussed in Section IV.6). However, it seems to be inconsistent with the observed spin flop transition in GdFe₂Si₂ [64], which implies that we may need to treat the correlation effect more precisely to understand the magnetisn in GdFe₂Si₂. On the other hand, GdRu₂Si₂ has a peak at $\bm{Q}\sim(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-$\bm{Q}$ states like a skyrmion lattice phase. For a spin model of GdRu₂Si₂, 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 $\bm{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 easy-axis 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₂Si₂. The value of $K$ in GdRu₂Si₂ is $\sim 0.43$ meV, which is in good agreement with the spin model simulation Note that the significance of weak easy-axis anisotropy can be straightforwardly understood: if the anisotropy is easy-plane type, a helical single-$\bm{Q}$ state is favored over multiple-$\bm{Q}$ states like a skyrmion phase, and if the easy-axis anisotropy is too strong, Ising-type magnetic order is favored over finite-$\bm{Q}$ structures. Thus, we conclude that the SDFT+$U$ result in GdRu₂Si₂ is fully consistent with the observation of the skyrmion lattice phase, and the present method correctly captures the essential physics behind the magnetic structure.

One of the advantages of evaluating $J(\bm{q})$ based on the tight-binding formalism is that it is easy to discuss the microscopic origin of the magnetic interaction [65, 48]. Although $J(\bm{q})$ was decomposed into the Gd-4$f$ and 5$d$ 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-4$f$ 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 interactions in this compound since it is included in the Gd-4$f$ 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-$\bm{Q}$ instability. Since the $\bm{q}$ dependence of $J(\bm{q})$ mainly comes from the Gd-5$d$ manifold according to Fig. 6(a), it would be interesting to decompose $J(\bm{q})$ into each 5$d$ orbital based on the SDFT+$U$ electronic structures.

Figure 6(b) and (c) show the diagonal and off-diagonal contributions to $J(\bm{q})$ in the Gd-5$d$ manifold, respectively. Among the five 5$d$ orbitals, the contributions from the $d_{z^{2}}$ and $d_{xy}$ orbitals are negligibly small, and the $\bm{q}$ dependence of $J(\bm{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 $\Gamma$ point, favoring the ferromagnetic ground state. In contrast, the other diagonal and off-diagonal ones show the peak at $\bm{Q}\sim(0.25,0,0)$ or almost flat along the $\Gamma$-X line. These results indicate that the orbital frustration exists in the 5$d$ orbital manifold, and whether the system favors the ferromagnetic or finite-$\bm{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₂Si₂ but also in GdFe₂Si₂ 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₂Si₂ than GdRu₂Si₂. This is mainly due to the suppression of the finite-$\bm{Q}$ peak in the contribution from the other $d$ orbitals in GdFe₂Si₂.

Here, let us look into the detail of the microscopic origin of $J(\bm{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-5$d$ orbitals, 4$f$ orbitals and their hybridization terms, respectively. Our calculations for Gd$T_{2}X_{2}$ show that, in the presence of $U$, the 5$d$ contribution always overcomes that of 4$f$ since the latter is strongly suppressed. This clearly indicates that both the super-exchange interaction $J_{\rm ex}^{f}\sim t_{f}^{2}/U_{f}$ and the conventional RKKY interaction $J_{\rm RKKY}\sim(V^{2}/U_{f})^{2}\chi^{d}$, where $\chi^{d}$ is the spin susceptibility in the 5$d$ orbitals, are negligibly small in these compounds. Conversely, the strong magnetic interaction from the 5$d$ 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 4$f$ orbitals. Here, we note that there are two possibilities for the origin of the exchange splitting in the 5$d$ states. One is the Coulomb interaction inherent in the 5$d$ orbitals, i.e., the second term in eq. (11), in which case the spin instability is described purely by the 5$d$ orbitals. The other is the magnetic coupling between the 5$d$ and 4$f$ 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 one-body calculation. Then, we can not distinguish these two contributions from the resulting exchange splittings. Investigating the microscopic origin of $J(\bm{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.

Let us move on to the results of $J(\bm{q})$ calculated with SDFT+$U$ for Gd$T_{2}$$X_{2}$ and Eu$T_{2}$$X_{2}$ in Fig. 9 and Fig. 9, respectively. From these figures, we see that the general trend of $J(\bm{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., 3$d$, 4$d$, or 5$d$) of the transition metals $T$ dramatically change the structure of $J({\bm{q}})$. For example, for the Gd$T_{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 Eu$T_{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. 9 and 9, the most promising candidates showing skyrmion lattice phase would be compounds with $T$ being the group 8 transition metal. In particular, GdRu₂Si₂, GdRu₂Ge₂, GdOs₂Si₂, GdOs₂Ge₂, and EuOs₂Si₂ have the highest peaks at finite $\bm{Q}$ along the $\Gamma$-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₂Si₂. 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 $\Gamma$-X line. However, the peaks in GdW${}_{2}X_{2}$ and EuCo${}_{2}X_{2}$ may be too small to stabilize the magnetic order.

Up to now, we have discussed the magnetic interaction $J(\bm{q})$ and the corresponding modulation vectors $\bm{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 $\bm{q}$-dependence far from the $\Gamma$ 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 $\bm{q}$-dependence of the total energies for Gd$T_{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(\bm{q})$ represents the most stable modulation vector, we can see that GdRu₂Si₂, GdRu₂Ge₂, and GdOs₂Si₂ favor the finite-$Q$ state with $\bm{Q}\sim(0.20,0,0)$-$(0.25,0,0)$. On the other hand, GdFe₂Si₂ and GdFe₂Ge₂ have the minimum at the $\Gamma$ 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₂Ge₂ does not show the finite-$\bm{Q}$ instability in contrast with the Liechtenstein result, the peak of $J(\bm{q})$ in GdOs₂Si₂ is more fragile than the others as is shown in Fig. 9(e). Thus, the disagreement with the Liechtenstein method is not serious. For the comparison, we also show $E(\bm{q})$ of Gd$T_{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-$\bm{Q}$ instability becomes stronger in all cases than in SDFT+$U$. Here, GdFe₂Si₂ shows the finite-$\bm{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(\bm{q})$ for the other candidates for showing the skyrmion phase predicted by the Liechtenstein calculations, namely, Gd$T_{2}X_{2}$ with $T=$ Re, Ag, Au and $X=$ Si, Ge, EuAg₂Si₂ and EuAg₂Ge₂. We can see that GdRe₂Si₂, and GdRe₂Ge₂ show nearly flat $\bm{q}$ dependence near the $\Gamma$ point, and thus, it is possible that some subtle perturbations select a finite-$\bm{Q}$ state. On the other hand, GdAg₂$X_{2}$, GdAu₂$X_{2}$, and EuAg₂$X_{2}$ have a broad peak at $\bm{Q}\sim(\frac{1}{2},0,0)$ corresponding to the 90-degree 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 $\bm{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.