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

Takuya Nomoto<EMAIL_ADDRESS>Research Center for Advanced Science
and Technology, University of Tokyo, Komaba, Meguro-ku, Tokyo 153-8904, Japan
Ryotaro Arita<EMAIL_ADDRESS>Research Center for Advanced Science and
Technology, University of Tokyo, Komaba, Meguro-ku, Tokyo 153-8904, Japan
RIKEN Center for Emergent Matter Science (CEMS), Wako 351-0198, Japan

###### Abstract

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(\bm{q})$ with a peak at finite-$\bm{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 Perspective, we
focus on recent developments of the first-principles evaluations of these
properties and apply them to the prototype systems Gd$T_{2}X_{2}$ and
Eu$T_{2}X_{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-${\bm{Q}}$
instability and show that competition among Gd-5$d$ orbitals determines
whether ferromagnetism or a finite-${\bm{Q}}$ structure is favored in
Gd$T_{2}$Si2 with $T=$ Fe and Ru.  Our systematic calculations reveal that
GdRu2$X_{2}$, GdOs2$X_{2}$, and GdRe${}_{2}X_{2}$ are promising, while
GdAg${}_{2}X_{2}$, GdAu${}_{2}X_{2}$, and EuAg${}_{2}X_{2}$ are possible
candidates as the skyrmion host materials. Analysis based on a spin spiral
calculation for the candidate materials is also presented.

††preprint: paper

##  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$ being
magnetic stiffness and $D$ being a DM constant. Once the spiral state is
stabilized, a triple-$\bm{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\ll 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 Gd2PdSi3 with
Gd forming a two-dimensional triangular lattice network [24, 25], followed by
Gd3Ru4Al12 with a Kagome lattice [26], GdRu2Si2 [27], EuAl4 [28, 29, 30, 31],
and EuGa4 [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 Mn1-xFexGe [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
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 momentum-space 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-$\bm{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
GdRu2Si2 [48, 49], Gd2PdSi3 [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-$\bm{Q}$ modulations [47]. Since these
effects on skyrmion formation are well-established, in this Perspective, we
aim to highlight recent progress in first-principles 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 Gd$T_{2}X_{2}$ and
Eu$T_{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
ThCr2Si2-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(\bm{q})$ based on the so-called Liechtenstein method with the
Wannier tight-binding model formalism. The calculated $J(\bm{q})$ is in good
agreement with $E(\bm{q})$ directly evaluated by the spin spiral calculations
within SDFT+$U$, supporting the validity of the present Liechtenstein
calculations. Based on the results, we argue that the finite-${\bm{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-5$d$ orbital contributions determines whether a ferromagnetism or
finite-${\bm{Q}}$ structure is favored in Gd$T_{2}$Si2 with $T=$ Fe and Ru.
According to our systematic calculations, GdRu2$X_{2}$, GdOs2$X_{2}$,
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 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(\bm{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,

$\displaystyle H=\sum_{12}(t_{12}+v_{12})c_{1}^{\dagger}c_{2},$ (1)

where indices 1 and 2 run over all degrees of freedom that specify the Wannier
functions, including atomic sites $i$, orbitals $\ell$, and spins $\sigma$.
$c_{1}^{\dagger}$ ($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,

$\displaystyle J_{ij}=-\frac{1}{2}{\rm
Tr}_{n\ell\sigma}[G_{ij}(\omega_{n})\Sigma^{i}G_{ji}(i\omega_{n})\Sigma^{j}].$
(2)

Here, ${\rm Tr}_{n\ell\sigma}=T\sum_{n}{\rm Tr}_{\ell\sigma}$, where ${\rm
Tr}_{\ell\sigma}$ denotes the trace over the orbital and spin spaces.
$\omega_{n}=\pi T(2n+1)$ denotes the Matsubara frequency, and $T$ is the
temperature. The Green’s function $G(\omega_{n})$ and the magnetic
perturbation matrix at $i$ site $\Sigma^{i}$ are defined by,

$\displaystyle[G(\omega_{n})]^{-1}_{12}$
$\displaystyle=i\omega_{n}\delta_{12}-t_{12}-v_{12}$ (3)
$\displaystyle\Sigma_{12}^{i}$
$\displaystyle=\tilde{v}_{i\ell_{1},i\ell_{2}}\sigma^{x}_{\sigma_{1}\sigma_{2}}$
(4)

where we approximate $\Sigma^{i}_{12}$ as a local quantity. $\sigma^{\mu}$
($\mu=x,y,z$) is the Pauli matrix, and
$\tilde{v}_{i_{1}\ell_{1},i_{2}\ell_{2}}$ is defined by
$v_{i_{1}\ell_{1}\sigma_{1},i_{2}\ell_{2}\sigma_{2}}=\tilde{v}_{i_{1}\ell_{1},i_{2}\ell_{i}}\sigma^{z}_{\sigma_{1}\sigma_{2}}$.
Then, $G_{ij}(\omega_{n})$ is defined as a sub-matrix of $G_{12}(\omega_{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(\bm{q})$ is obtained as follows:

$\displaystyle
J(\bm{q})=\frac{1}{N}\sum_{ij}e^{-i\bm{q}\cdot(\bm{R}_{i}-\bm{R}_{j})}J_{ij}-J_{0},$
(5)

where $\bm{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 Gd$T_{2}X_{2}$ and Eu$T_{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(\bm{q})$,
which is essential for the SDFT+$U$ calculations. Based on the obtained
$J(\bm{q})$, we can see whether the system favors the finite-$\bm{Q}$
modulation.

On the other hand, we can evaluate the energy of the spin spiral states with
the modulation vector $\bm{Q}$, denoted by $E(\bm{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(\bm{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-$\bm{Q}$ structure, which makes the calculation with
small $\bm{Q}$ tractable. However, despite these advantages, it requires a
calculation for every $\bm{Q}$, and thus, is much more expensive than the
Liechtenstein method if one wants to see the detailed structures of
$E(\bm{q})$. Thus, we use these two methods in a complementary manner.

##  III Calculation details

For the evaluations of $J(\bm{q})$ and $E(\bm{q})$, we use Vienna ab initio
simulation package (VASP) [59] for the electronic structure calculations.
Here, the exchange-correlation 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 Gd$T_{2}X_{2}$ and
Eu$T_{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\times 16\times 16$ Monkhorst-Pack
$\bm{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\times 21\times 21$ Monkhorst-Pack
$\bm{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\times 16\times 16$
Monkhorst-Pack $\bm{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\times 10\times 10$ uniform $\bm{k}$-grid 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(\bm{q})$ are evaluated by Eqs. (2) and (5), respectively. Here, we set
the temperature $T=58$ K and employ $48\times 48\times 48$ uniform
$\bm{k}$-grid.

##  IV Results and Discussion

###  IV.1 Magnetocrystalline anisotropy

Figure 1:  Magnetocrystalline anisotropy of GdRu2Si2 and GdRu2Ge2. The label
“without $U$” stands for the result of the standard SDFT calculation. The
others correspond to SDFT+$U$ calculations with $U=6.7$ eV.

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 Gd3+ and Eu2+, 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 Gd3+\- and Eu2+-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.

Figure 2:  Magnetocrystalline anisotropy of (a) Gd$T_{2}$Si2, (b)
Gd$T_{2}$Ge2, (c) Eu$T_{2}$Si2 and (d) Eu$T_{2}$Ge2. 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.

As an example, let us first see in Fig. 1 the calculated magnetocrystalline
anisotropy,

$\displaystyle\Delta E:=E(\bm{m}//\bm{x})-E(\bm{m}//\bm{z})$ (6)

for the skyrmion compound GdRu2Si2 and a reference compound GdRu2Ge2 [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 Gd3+, 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 GdRu2Si2 and
GdRu2Ge2. 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 ThCr2Si2-type
crystal structure.

###  IV.2 Fermi surface and magnetic instability

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 Gd3+ and Eu2+ 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 EuMn2Si2, EuZr2Si2,
EuZr2Ge2, EuTc2Si2, EuTc2Ge2, EuHf2Ge2, EuTa2Si2 and EuRe2Ge2. 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.

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

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 GdRu2Si2 and GdFe2Si2 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 Gd3+ and Eu2+.

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,

$\displaystyle\chi_{nn^{\prime}}({\bm{q}})=-\frac{1}{N}\sum_{\bm{k}}\frac{f(\varepsilon_{n{\bm{k}}+{\bm{q}}})-f(\varepsilon_{n^{\prime}{\bm{k}}})}{\varepsilon_{n{\bm{k}}+{\bm{q}}}-\varepsilon_{n^{\prime}{\bm{k}}}},$
(7)

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,

$\displaystyle\chi_{\mu\nu}({\bm{q}})$
$\displaystyle=\frac{1}{N}\sum_{\bm{R}}\sum_{1234}e^{i\phi_{13}^{\bm{q}}(\bm{R})}s^{\mu}_{13}s^{\nu}_{24}\chi_{13,24}(\bm{R}),$
(8) $\displaystyle\chi_{13,24}(\bm{R})$
$\displaystyle=-T\sum_{\omega_{n}}G_{12}(i\omega_{n},\bm{R})G_{43}(i\omega_{n},-\bm{R}).$
(9)

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.

###  IV.3 $J(\bm{q})$ in GdRu2Si2 and GdFe2Si2 evaluated by the Liechtenstein
method

Figure 4:  Modulation vector $\bm{Q}$ dependence of the spin exchange
interaction $J(\bm{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.

Figures 4(a) and (b) show $J(\bm{q})$ for GdRu2Si2 and GdFe2Si2 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 GdRu2Si2 favors the spin spiral state with $\bm{Q}\sim(0.25,0,0)$
both in SDFT and SDFT+$U$. On the other hand, in GdFe2Si2, 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
GdFe2Si2 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 GdFe2Si2 [64], which implies that we may need to treat the
correlation effect more precisely to understand the magnetisn in GdFe2Si2. On
the other hand, GdRu2Si2 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 GdRu2Si2, an
analysis based on the Landau-Lifshitz-Gilbert (LLG) equation was performed in
Ref. [49], using the following model for the internal energy,

$\displaystyle
E=-\frac{1}{2}\sum_{ij}J_{ij}{\bm{m}}_{i}\cdot{\bm{m}}_{j}-K\sum_{i}(\bm{m}_{i}\cdot\bm{e}_{z})^{2}-\bm{B}\cdot\sum_{i}\bm{m}_{i},$
(10)

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 GdRu2Si2. The
value of $K$ in GdRu2Si2 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 GdRu2Si2 is fully consistent with the
observation of the skyrmion lattice phase, and the present method correctly
captures the essential physics behind the magnetic structure.

Figure 5:  (a) A magnetic phase diagram of GdRu2Si2 in the $(K,B)$ parameter
space, where $K$ is the magnetocrystalline anisotropy and $B$ is the magnetic
field. (b)-(d) Representative magnetic textures described by the Gd moments.
(b) Helical state for $(K,B)=(0.2,0.6)$, (c) Skyrmion state for
$(K,B)=(0.2,1.0)$, and (d) cycloidal state for $(K,B)=(0.1,1.0)$. Reproduced
from Ref. [49].

###  IV.4 Orbital decomposition of $J(\bm{q})$ calculated based on the
Liechtenstein method

Figure 6:  Orbital decomposed $J(\bm{q})$ calculated by the Liechtenstein
method. (a) The contributions from the Gd-5$d$ and 4$f$ orbitals. (b) The
diagonal contribution to in the Gd-5$d$ manifold. Here, the indices 1 to 5
correspond to the $d_{z^{2}},d_{xz},d_{yz},d_{x^{2}-y^{2}}$, and $d_{xy}$
orbitals, respectively. (c) The off-diagonal contributions in the Gd-5$d$
manifold.

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 7:  Orbital decomposed $J(\bm{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.

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 GdRu2Si2 but also in GdFe2Si2 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 GdFe2Si2 than GdRu2Si2. This is mainly due to the suppression of the
finite-$\bm{Q}$ peak in the contribution from the other $d$ orbitals in
GdFe2Si2.

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,

$\displaystyle H_{d}$
$\displaystyle=t_{d}\sum_{\braket{i,j}\sigma}d^{\dagger}_{i\sigma}d_{j\sigma}+U_{d}\sum_{i}n_{i\uparrow}^{d}n_{i\downarrow}^{d},$
(11) $\displaystyle H_{f}$
$\displaystyle=t_{f}\sum_{\braket{i,j}\sigma}f^{\dagger}_{i\sigma}f_{j\sigma}+U_{f}\sum_{i}n_{i\uparrow}^{f}n_{i\downarrow}^{f},$
(12) $\displaystyle H_{df}$
$\displaystyle=V\sum_{i\sigma}(d^{\dagger}_{i\sigma}f_{i\sigma}+f^{\dagger}_{i\sigma}d_{i\sigma})+J_{cf}\sum_{i}{\bm{s}}^{d}\cdot{\bm{s}}^{f}.$
(13)

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.

Figure 8:  Spin interactions $J(\bm{q})$ in Gd$T_{2}$$X_{2}$ 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
3$d$, 4$d$, and 5$d$ elements as $T$, respectively. The solid (dashed) lines
with filled (open) symbols stand for $X=$ Si (Ge). Note that the high symmetry
points Z, $\Gamma$, and X correspond to $(0,0,\frac{2\pi}{c}),(0,0,0)$ and
$(\frac{2\pi}{a},0,0)$ in the cartesian frame of the reciprocal lattice space,
respectively.

Figure 9:  Spin interactions $J(\bm{q})$ in Eu$T_{2}$$X_{2}$ calculated by the
Liechtenstein method based on the SDFT+$U$ electronic structures. Color, line
and point types are the same as Fig. 9.

###  IV.5 Systematic evaluation of $J(\bm{q})$ based on the Liechtenstein
method

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, GdRu2Si2, GdRu2Ge2, GdOs2Si2, GdOs2Ge2, and EuOs2Si2 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 GdRu2Si2. 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.

###  IV.6 Modulation vector based on spin spiral calculations

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:  Modulation vector $\bm{q}$ dependence of the ground state energy
$E(\bm{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 Gd$T_{2}X_{2}$ with
$T=$ Fe, Ru, Os and $X=$ Si, Ge. (b) The results based on the SDFT for
Gd$T_{2}X_{2}$ 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 Gd$T_{2}X_{2}$ with $T=$ Re,
Ag, Au and $X=$ Si, Ge.

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 GdRu2Si2, GdRu2Ge2, and GdOs2Si2 favor the finite-$Q$
state with $\bm{Q}\sim(0.20,0,0)$-$(0.25,0,0)$. On the other hand, GdFe2Si2
and GdFe2Ge2 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 GdOs2Ge2
does not show the finite-$\bm{Q}$ instability in contrast with the
Liechtenstein result, the peak of $J(\bm{q})$ in GdOs2Si2 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, GdFe2Si2 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, EuAg2Si2 and EuAg2Ge2. We can see that GdRe2Si2, and GdRe2Ge2
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, GdAg2$X_{2}$, GdAu2$X_{2}$, and EuAg2$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.

##  V Conclusion and outlook

In this Perspective, we present systematic first-principles calculations for
Gd$T_{2}X_{2}$ and Eu$T_{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(\bm{q})$ based on the Liechtenstein method and show that the
obtained $J(\bm{q})$ is in good agreement with $E(\bm{q})$ by the spin spiral
calculations. Our calculations indicate that the finite-${\bm{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-${\bm{Q}}$ structure is favored in Gd$T_{2}$Si2 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-$\bm{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.

##  VI Acknowledgments

This work was supported by JSPS-KAKENHI (No. 22H00290, 21H04437, 21H04990, and
19H05825) and JST-PRESTO (No. JPMJPR20L7). We also acknowledge the Center for
Computational Materials Science, Institute for Materials Research, Tohoku
University, for the use of MASAMUNE-IMR (Project No. 202112-SCKXX-0010).

##  VII AUTHOR DECLARATIONS

Conflict of Interest

The authors have no conflicts to disclose.

Author Contributions

Takuya Nomoto:Formal analysis (lead); Writing - original draft. Ryotaro Arita:
Formal analysis (supporting); Writing - review & editing.

DATA AVAILABILITY

The data that support the findings of this study are available from the
corresponding author upon reasonable request.

## References

  * Nagaosa and Tokura [2013] N. Nagaosa and Y. Tokura, Topological properties and dynamics of magnetic skyrmions, Nature Nanotechnology 8, 899 (2013).
  * Göbel  _et al._ [2021] B. Göbel, I. Mertig, and O. A. Tretiakov, Beyond skyrmions: Review and perspectives of alternative magnetic quasiparticles, Physics Reports 895, 1 (2021), beyond skyrmions: Review and perspectives of alternative magnetic quasiparticles.
  * Tokura and Kanazawa [2021] Y. Tokura and N. Kanazawa, Magnetic skyrmion materials, Chemical Reviews 121, 2857 (2021), pMID: 33164494.
  * Finocchio  _et al._ [2016] G. Finocchio, F. B. ttner, R. Tomasello, M. Carpentieri, and M. Kläui, Magnetic skyrmions: from fundamental to applications, Journal of Physics D: Applied Physics 49, 423001 (2016).
  * Fert  _et al._ [2017] A. Fert, N. Reyren, and V. Cros, Magnetic skyrmions: advances in physics and potential applications, Nature Reviews Materials 2, 17031 (2017).
  * Luo and You [2021] S. Luo and L. You, Skyrmion devices for memory and logic applications, APL Materials 9, 050901 (2021).
  * Marrows and Zeissler [2021] C. H. Marrows and K. Zeissler, Perspective on skyrmion spintronics, Applied Physics Letters 119, 250502 (2021).
  * Wang  _et al._ [2022] H. Wang, Y. Dai, G.-M. Chow, and J. Chen, Topological hall transport: Materials, mechanisms and potential applications, Progress in Materials Science 130, 100971 (2022).
  * Parkin  _et al._ [2008] S. S. P. Parkin, M. Hayashi, and L. Thomas, Magnetic domain-wall racetrack memory, Science 320, 190 (2008).
  * Sampaio  _et al._ [2013] J. Sampaio, V. Cros, S. Rohart, A. Thiaville, and A. Fert, Nucleation, stability and current-induced motion of isolated magnetic skyrmions in nanostructures, Nature Nanotechnology 8, 839 (2013).
  * Zhang  _et al._ [2018] X. Zhang, W. Cai, X. Zhang, Z. Wang, Z. Li, Y. Zhang, K. Cao, N. Lei, W. Kang, Y. Zhang, H. Yu, Y. Zhou, and W. Zhao, Skyrmions in magnetic tunnel junctions, ACS Applied Materials & Interfaces 10, 16887 (2018).
  * Huang  _et al._ [2017] Y. Huang, W. Kang, X. Zhang, Y. Zhou, and W. Zhao, Magnetic skyrmion-based synaptic devices, Nanotechnology 28, 08LT02 (2017).
  * Song  _et al._ [2020] K. M. Song, J.-S. Jeong, B. Pan, X. Zhang, J. Xia, S. Cha, T.-E. Park, K. Kim, S. Finizio, J. Raabe, J. Chang, Y. Zhou, W. Zhao, W. Kang, H. Ju, and S. Woo, Skyrmion-based artificial synapses for neuromorphic computing, Nature Electronics 3, 148 (2020).
  * Lin  _et al._ [1973] Y. S. Lin, P. J. Grundy, and E. A. Giess, Bubble domains in magnetostatically coupled garnet films, Applied Physics Letters 23, 485 (1973).
  * Garel and Doniach [1982] T. Garel and S. Doniach, Phase transitions with spontaneous modulation-the dipolar ising ferromagnet, Phys. Rev. B 26, 325 (1982).
  * Mühlbauer  _et al._ [2009] S. Mühlbauer, B. Binz, F. Jonietz, C. Pfleiderer, A. Rosch, A. Neubauer, R. Georgii, and P. Böni, Skyrmion lattice in a chiral magnet, Science 323, 915 (2009).
  * Bogdanov and Yablonskii [1989] A. N. Bogdanov and D. A. Yablonskii, Thermodynamically stable “vortices” in magnetically ordered crystals. the mixed state of magnets, Sov. Phys. JETP 68, 101 (1989).
  * Rößler  _et al._ [2006] U. K. Rößler, A. N. Bogdanov, and C. Pfleiderer, Spontaneous skyrmion ground states in magnetic metals, Nature 442, 797 (2006).
  * Fert and Levy [1980] A. Fert and P. M. Levy, Role of anisotropic exchange interactions in determining the properties of spin-glasses, Phys. Rev. Lett. 44, 1538 (1980).
  * Heide  _et al._ [2008] M. Heide, G. Bihlmayer, and S. Blügel, Dzyaloshinskii-moriya interaction accounting for the orientation of magnetic domains in ultrathin films: Fe/w(110), Phys. Rev. B 78, 140403 (2008).
  * Freimuth  _et al._ [2014] F. Freimuth, S. B. gel, and Y. Mokrousov, Berry phase theory of dzyaloshinskii– moriya interaction and spin– orbit torques, Journal of Physics: Condensed Matter 26, 104202 (2014).
  * Kikuchi  _et al._ [2016] T. Kikuchi, T. Koretsune, R. Arita, and G. Tatara, Dzyaloshinskii-moriya interaction as a consequence of a doppler shift due to spin-orbit-induced intrinsic spin current, Phys. Rev. Lett. 116, 247201 (2016).
  * Koretsune  _et al._ [2018] T. Koretsune, T. Kikuchi, and R. Arita, First-principles evaluation of the dzyaloshinskii– moriya interaction, Journal of the Physical Society of Japan 87, 041011 (2018).
  * Saha  _et al._ [1999] S. R. Saha, H. Sugawara, T. D. Matsuda, H. Sato, R. Mallik, and E. V. Sampathkumaran, Magnetic anisotropy, first-order-like metamagnetic transitions, and large negative magnetoresistance in single-crystal ${\mathrm{gd}}_{2}{\mathrm{pdsi}}_{3}$, Phys. Rev. B 60, 12162 (1999).
  * Kurumaji  _et al._ [2019] T. Kurumaji, T. Nakajima, M. Hirschberger, A. Kikkawa, Y. Yamasaki, H. Sagayama, H. Nakao, Y. Taguchi, T. hisa Arima, and Y. Tokura, Skyrmion lattice with a giant topological hall effect in a frustrated triangular-lattice magnet, Science 365, 914 (2019).
  * Hirschberger  _et al._ [2019] M. Hirschberger, T. Nakajima, S. Gao, L. Peng, A. Kikkawa, T. Kurumaji, M. Kriener, Y. Yamasaki, H. Sagayama, H. Nakao, K. Ohishi, K. Kakurai, Y. Taguchi, X. Yu, T.-h. Arima, and Y. Tokura, Skyrmion phase and competing magnetic orders on a breathing kagomé lattice, Nature Communications 10, 5831 (2019).
  * Khanh  _et al._ [2020] N. D. Khanh, T. Nakajima, X. Yu, S. Gao, K. Shibata, M. Hirschberger, Y. Yamasaki, H. Sagayama, H. Nakao, L. Peng, K. Nakajima, R. Takagi, T.-h. Arima, Y. Tokura, and S. Seki, Nanometric square skyrmion lattice in a centrosymmetric tetragonal magnet, Nature Nanotechnology 15, 444 (2020).
  * Shang  _et al._ [2021] T. Shang, Y. Xu, D. J. Gawryluk, J. Z. Ma, T. Shiroka, M. Shi, and E. Pomjakushina, Anomalous hall resistivity and possible topological hall effect in the ${\mathrm{eual}}_{4}$ antiferromagnet, Phys. Rev. B 103, L020405 (2021).
  * Kaneko  _et al._ [2021] K. Kaneko, T. Kawasaki, A. Nakamura, K. Munakata, A. Nakao, T. Hanashima, R. Kiyanagi, T. Ohhara, M. Hedo, T. Nakama, and Y. Ōnuki, Charge-density-wave order and multiple magnetic transitions in divalent europium compound eual4, Journal of the Physical Society of Japan 90, 064704 (2021).
  * Zhu  _et al._ [2022] X. Y. Zhu, H. Zhang, D. J. Gawryluk, Z. X. Zhen, B. C. Yu, S. L. Ju, W. Xie, D. M. Jiang, W. J. Cheng, Y. Xu, M. Shi, E. Pomjakushina, Q. F. Zhan, T. Shiroka, and T. Shang, Spin order and fluctuations in the ${\mathrm{eual}}_{4}$ and ${\mathrm{euga}}_{4}$ topological antiferromagnets: A $\mu\mathrm{SR}$ study, Phys. Rev. B 105, 014423 (2022).
  * Takagi  _et al._ [2022] R. Takagi, N. Matsuyama, V. Ukleev, L. Yu, J. S. White, S. Francoual, J. R. L. Mardegan, S. Hayami, H. Saito, K. Kaneko, K. Ohishi, Y. Ōnuki, T.-h. Arima, Y. Tokura, T. Nakajima, and S. Seki, Square and rhombic lattices of magnetic skyrmions in a centrosymmetric binary compound, Nature Communications 13, 1472 (2022).
  * Zhang  _et al._ [2021] H. Zhang, X. Y. Zhu, Y. Xu, D. J. Gawryluk, W. Xie, S. L. Ju, M. Shi, T. Shiroka, Q. F. Zhan, E. Pomjakushina, and T. Shang, Giant magnetoresistance and topological hall effect in the euga4 antiferromagnet, Journal of Physics: Condensed Matter 34, 034005 (2021).
  * [33] Y. Ōnuki, M. Kakihana, W. Iha, K. Nakaima, D. Aoki, A. Nakamura, F. Honda, M. Nakashima, Y. Amako, J. Gouchi, Y. Uwatoko, S. Nakamura, T. Sakakibara, T. Takeuchi, Y. Haga, H. Ikeda, H. Harima, M. Hedo, and T. Nakama, Unique skyrmion phases and conduction electrons in cubic chiral antiferromagnet euptsi and related compounds, in  _Proceedings of the International Conference on Strongly Correlated Electron Systems (SCES2019)_.
  * Sakakibara  _et al._ [2021] T. Sakakibara, S. Nakamura, S. Kittaka, M. Kakihana, M. Hedo, T. Nakama, and Y. Ōnuki, Magnetic phase transitions of the 4f skyrmion compound euptsi studied by magnetization measurements, Journal of the Physical Society of Japan 90, 064701 (2021).
  * Kanazawa  _et al._ [2011] N. Kanazawa, Y. Onose, T. Arima, D. Okuyama, K. Ohoyama, S. Wakimoto, K. Kakurai, S. Ishiwata, and Y. Tokura, Large topological hall effect in a short-period helimagnet mnge, Phys. Rev. Lett. 106, 156603 (2011).
  * Shibata  _et al._ [2013] K. Shibata, X. Z. Yu, T. Hara, D. Morikawa, N. Kanazawa, K. Kimoto, S. Ishiwata, Y. Matsui, and Y. Tokura, Towards control of the size and helicity of skyrmions in helimagnetic alloys by spin– orbit coupling, Nature Nanotechnology 8, 723 (2013).
  * Hamamoto  _et al._ [2015] K. Hamamoto, M. Ezawa, and N. Nagaosa, Quantized topological hall effect in skyrmion crystal, Phys. Rev. B 92, 115417 (2015).
  * Matsui  _et al._ [2021] A. Matsui, T. Nomoto, and R. Arita, Skyrmion-size dependence of the topological hall effect: A real-space calculation, Phys. Rev. B 104, 174432 (2021).
  * Okubo  _et al._ [2012] T. Okubo, S. Chung, and H. Kawamura, Multiple-$q$ states and the skyrmion lattice of the triangular-lattice heisenberg antiferromagnet under magnetic fields, Phys. Rev. Lett. 108, 017206 (2012).
  * Leonov and Mostovoy [2015] A. O. Leonov and M. Mostovoy, Multiply periodic states and isolated skyrmions in an anisotropic frustrated magnet, Nature Communications 6, 8275 (2015).
  * Chen  _et al._ [2016] J. P. Chen, D.-W. Zhang, and J.-M. Liu, Exotic skyrmion crystals in chiral magnets with compass anisotropy, Scientific Reports 6, 29126 (2016).
  * Chen  _et al._ [2018] J. Chen, D.-W. Zhang, Y. Chen, X. Gao, and J.-M. Liu, Compass-anisotropy-modulated helical states and skyrmion crystals in chiral magnets, Physics Letters A 382, 2944 (2018).
  * Heinze  _et al._ [2011] S. Heinze, K. von Bergmann, M. Menzel, J. Brede, A. Kubetzka, R. Wiesendanger, G. Bihlmayer, and S. Blü gel, Spontaneous atomic-scale magnetic skyrmion lattice in two dimensions, Nature Physics 7, 713 (2011).
  * Hayami  _et al._ [2017] S. Hayami, R. Ozawa, and Y. Motome, Effective bilinear-biquadratic model for noncoplanar ordering in itinerant magnets, Phys. Rev. B 95, 224424 (2017).
  * Ozawa  _et al._ [2016] R. Ozawa, S. Hayami, K. Barros, G.-W. Chern, Y. Motome, and C. D. Batista, Vortex crystals with chiral stripes in itinerant magnets, Journal of the Physical Society of Japan 85, 103703 (2016).
  * Ozawa  _et al._ [2017] R. Ozawa, S. Hayami, and Y. Motome, Zero-field skyrmions with a high topological number in itinerant magnets, Phys. Rev. Lett. 118, 147205 (2017).
  * Hayami and Motome [2021] S. Hayami and Y. Motome, Topological spin crystals by itinerant frustration, Journal of Physics: Condensed Matter 33, 443001 (2021).
  * Nomoto  _et al._ [2020a] T. Nomoto, T. Koretsune, and R. Arita, Formation mechanism of the helical $\bm{q}$ structure in gd-based skyrmion materials, Phys. Rev. Lett. 125, 117204 (2020a).
  * Bouaziz  _et al._ [2022] J. Bouaziz, E. Mendive-Tapia, S. Blügel, and J. B. Staunton, Fermi-surface origin of skyrmion lattices in centrosymmetric rare-earth intermetallics, Phys. Rev. Lett. 128, 157206 (2022).
  * Grytsiuk  _et al._ [2020] S. Grytsiuk, J.-P. Hanke, M. Hoffmann, J. Bouaziz, O. Gomonay, G. Bihlmayer, S. Lounis, Y. Mokrousov, and S. Blü gel, Topological– chiral magnetic interactions driven by emergent orbital magnetism, Nature Communications 11, 511 (2020).
  * Mendive-Tapia  _et al._ [2021] E. Mendive-Tapia, M. dos Santos Dias, S. Grytsiuk, J. B. Staunton, S. Blügel, and S. Lounis, Short period magnetization texture of b20-mnge explained by thermally fluctuating local moments, Phys. Rev. B 103, 024410 (2021).
  * Korotin  _et al._ [2015] D. M. Korotin, V. V. Mazurenko, V. I. Anisimov, and S. V. Streltsov, Calculation of exchange constants of the heisenberg model in plane-wave-based methods using the green’s function approach, Phys. Rev. B 91, 224405 (2015).
  * Nomoto  _et al._ [2020b] T. Nomoto, T. Koretsune, and R. Arita, Local force method for the ab initio tight-binding model: Effect of spin-dependent hopping on exchange interactions, Phys. Rev. B 102, 014444 (2020b).
  * Mostofi  _et al._ [2008] A. A. Mostofi, J. R. Yates, Y.-S. Lee, I. Souza, D. Vanderbilt, and N. Marzari, wannier90: A tool for obtaining maximally-localised wannier functions, Computer Physics Communications 178, 685 (2008).
  * Marzari  _et al._ [2012] N. Marzari, A. A. Mostofi, J. R. Yates, I. Souza, and D. Vanderbilt, Maximally localized wannier functions: Theory and applications, Rev. Mod. Phys. 84, 1419 (2012).
  * Chikano  _et al._ [2019] N. Chikano, K. Yoshimi, J. Otsuki, and H. Shinaoka, irbasis: Open-source database and software for intermediate-representation basis functions of imaginary-time green’ s function, Computer Physics Communications 240, 181 (2019).
  * Li  _et al._ [2020] J. Li, M. Wallerberger, N. Chikano, C.-N. Yeh, E. Gull, and H. Shinaoka, Sparse sampling approach to efficient ab initio calculations at finite temperature, Phys. Rev. B 101, 035144 (2020).
  * Sandratskii [1998] L. M. Sandratskii, Noncollinear magnetism in itinerant-electron systems: Theory and applications, Advances in Physics 47, 91 (1998).
  * Kresse and Furthmüller [1996] G. Kresse and J. Furthmüller, Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set, Phys. Rev. B 54, 11169 (1996).
  * Perdew  _et al._ [1996] J. P. Perdew, K. Burke, and M. Ernzerhof, Generalized gradient approximation made simple, Phys. Rev. Lett. 77, 3865 (1996).
  * Blöchl [1994] P. E. Blöchl, Projector augmented-wave method, Phys. Rev. B 50, 17953 (1994).
  * Kresse and Joubert [1999] G. Kresse and D. Joubert, From ultrasoft pseudopotentials to the projector augmented-wave method, Phys. Rev. B 59, 1758 (1999).
  * Garnier  _et al._ [1995] A. Garnier, D. Gignoux, N. Iwata, D. Schmitt, T. Shigeoka, and F. Zhang, Anisotropic metamagnetism in gdru2si2, Journal of Magnetism and Magnetic Materials 140-144, 899 (1995), international Conference on Magnetism.
  * Mihalik  _et al._ [2005] M. Mihalik, P. Svoboda, M. Mihalik, J. Vejpravová, and V. Sechovský, Magnetic and transport properties of gdfe2si2 single-crystal, Journal of Magnetism and Magnetic Materials 290-291, 606 (2005), proceedings of the Joint European Magnetic Symposia (JEMS’ 04).
  * Kvashnin  _et al._ [2016] Y. O. Kvashnin, R. Cardias, A. Szilva, I. Di Marco, M. I. Katsnelson, A. I. Lichtenstein, L. Nordström, A. B. Klautau, and O. Eriksson, Microscopic origin of heisenberg and non-heisenberg exchange interactions in ferromagnetic bcc fe, Phys. Rev. Lett. 116, 217202 (2016).
  * [66] J. Bouaziz and S. Blügel, in private communication.
  * Katsnelson and Lichtenstein [2000] M. I. Katsnelson and A. I. Lichtenstein, First-principles calculations of magnetic interactions in correlated systems, Phys. Rev. B 61, 8906 (2000).