Skyrmion crystal in the RKKY system on the two-dimensional triangular lattice

Kota Mitsumoto Hikaru Kawamura Molecular Photoscience Research Center, Kobe University, Kobe 657-8501, Japan

Abstract

We study the ordering properties of the isotropic RKKY Heisenberg model on the two-dimensional (2D) triangular lattice by extensive Monte Carlo simulations to get insights into the chiral-degenerate skyrmion crystal (SkX) of metallic magnets. Our Hamiltonian contains only the spin-quadratic RKKY interaction derived from the spherical Fermi surface, containing neither the nesting nor the many-body interaction. The SkX phase is stabilized under applied fields where the frustration associated with the oscillating nature of the RKKY interaction and the emergent many-body interactions generated by thermal fluctuations play important roles. Replica symmetry breaking, reported in our recent study on the 3D RKKY model [Phys. Rev. B 104, 184432 (2021)], turns out to be absent in the present 2D model. Implications to the SkX formation mechanism are discussed.

I Introduction

Topologically protected spin textures in magnets have attracted much recent attention. In particular, skyrmion, a swirling noncoplanar spin texture whose constituent spin directions wrap a sphere in spin space, has been studied quite extensively. Skyrmion is usually stabilized as a periodic array called the skyrmion crystal (SkX). The SkX state was first observed in chiral ferromagnets such as MnSi Mühlbauer et al. (2009); Neubauer et al. (2009), FeCoSi Münzer et al. (2010); Yu et al. (2010), and FeGe Yu et al. (2011) under magnetic fields, where the SkX state is induced by the anti-symmetric Dzyaloshinskii-Moriya (DM) interaction. A recent study by Okubo, Chung, and Kawamura Okubo et al. (2012) of the $J_{1}$ - $J_{3}$ ( $J_{2}$ ) Heisenberg model on the two-dimensional (2D) triangular lattice has revealed that the SkX state can also be realized in centrosymmetric magnets without the DM interaction. This model shows the single- $q$ incommensurate spiral phase in zero field and exhibits the triple- $q$ SkX phase under applied fields at finite temperatures. The SkX state is induced by frustration among nearest- and further-neighbor exchange interactions and stabilized by magnetic fields and thermal fluctuations. Nowadays, the centrosymmetric SkX systems have attracted much interest both experimentally Kurumaji et al. (2019); Takahashi et al. (2020) and theoretically Leonov and Mostovoy (2015); Hayami et al. (2016); Ozawa et al. (2017); Hayami et al. (2017); Lin and Batista (2018); Hayami and Motome (2019); Wang et al. (2020); Hayami and Yambe (2020); Hayami and Motome (2021a, b); Yambe and Hayami (2021); Wang et al. (2021).

The skyrmion is characterized by the integer topological charge in units of the solid angle $4\pi$ , whose sign represents the swirling direction of the skyrmion, which is also represented by the sign of the scalar spin chirality. When the DM interaction induces the SkX, the total chirality is determined to be a definite sign depending on materials. In contrast, in the frustration-induced SkX in centrosymmetric magnets, the total chirality can take both positive and negative values, i.e., either SkX or anti-SkX being selected by the spontaneous symmetry breaking. We call such a symmetric state a chiral-degenerate SkX state. The chiral-degenerate SkX might exhibit an interesting electromagnetic response, e.g., the topological Hall effect of both signs. Such a chiral-degenerate nature might also give rise to a new interesting phase called the $Z$ phase, a random domain phase consisting of both SkX and anti-SkX.

Some candidate materials of the frustration-induced SkX are metallic compounds, Gd₂PdSi₃ Saha et al. (1999); Kurumaji et al. (2019) and EuCuSb Takahashi et al. (2020). The primary interaction between localized magnetic moments is the Ruderman-Kittel-Kasuya-Yosida (RKKY) interaction mediated by conduction electrons Ruderman and Kittel (1954); Kasuya (1956); Yosida (1957). When the exchange coupling $J_{sd}$ between the conduction electron and the localized spin is sufficiently weak compared with the Fermi energy $\epsilon_{\rm F}$ , the RKKY interaction can be derived via the second-order perturbation with respect to $J_{sd}/\epsilon_{\rm F}$ . The RKKY interaction in metals oscillates in sign with the distance $r$ , providing frustration just as in the competing $J_{1}$ - $J_{3}$ ( $J_{1}$ - $J_{2}$ ) interaction in insulators, while it is the long-range interaction falling as $1/r^{3}$ in contrast to the short-range interaction. The similarity in the inherent frustration of the interaction suggests that the chiral-degenerate SkX might also be stabilized in RKKY metals.

The theoretical studies on metallic systems, most of which dealt with the $T=0$ case, have suggested that the SkX cannot be realized in centrosymmetric metals by only the isotropic RKKY interaction. Refs. Ozawa et al. (2017); Hayami et al. (2017); Hayami and Yambe (2020) investigated the SkX phase by numerical simulation on the Kondo lattice model and the analysis of the effective classical spin model in $q$ -space, and suggested that incommensurate wavenumbers dictated by the Fermi surface, i.e., nesting, and the four-body spin interactions arising from the higher-order perturbation, especially the biquadratic interaction with the positive coefficient, are essential in stabilizing the SkX. Ref. Wang et al. (2020) investigated the RKKY system mediated by 2D electron gas by the $T=0$ variational method, showing that the SkX state needed an easy-axis anisotropy to be stabilized, in addition to the frustration of the RKKY interaction. Meanwhile, these studies mostly dealt with the $T=0$ case, and the studies of the possible SkX state in metals have been rather scarce at finite temperatures. Since the SkX phase in the isotropic short-range system is stabilized by thermal fluctuations, it is important to explore the possible SkX phase of the RKKY system at finite temperatures.

In the present paper, we investigate by extensive Monte Carlo (MC) simulations the ordering properties of the isotropic classical Heisenberg model on the 2D triangular lattice interacting via the standard RKKY interaction associated with the spherical Fermi surface. Then, our model does not possess the nesting, the four-body (biquadratic) interaction, or the magnetic anisotropy. Our first aim is to examine whether the SkX state is ever possible, even in such a simple model of metallic systems.

Our very recent study on the 3D RKKY system Mitsumoto and Kawamura (2021) has indicated that the SkX is stabilized, but the state accompanies quite an exotic property of the replica-symmetry breaking (RSB) familiar in glassy systems Mézard et al. (1987); Mydosh (1993); Kawamura and Taniguchi (2015) but rarely seen in regularly ordered states Mitsumoto and Kawamura (2021). Hence, the SkX state in the 3D RKKY system is quite different in nature from the chiral-degenerate SkX phase observed in the 2D short-range model Okubo et al. (2012). The second aim of the present paper is to clarify whether an RSB phenomenon occurs in the 2D RKKY system.

By extensive MC simulations, we have shown that the 2D RKKY system certainly exhibits a chiral-degenerate triple- $q$ SkX phase under magnetic fields, which, however, does not accompany the RSB. The results demonstrate that neither the nesting, the four-body interaction, nor the magnetic anisotropy is indispensable for the SkX formation. Implications of the results to the SkX formation mechanism will then be discussed.

The rest of the paper is organized as follows. In Sec. II, we introduce our model and explain the MC simulation method employed. We present the results of our MC simulations in Sec. III, including the temperature versus magnetic-field phase diagram in Sec. III.1, spin structure factors and real-space spin and chirality configurations in Sec. III.2, temperature dependences of several physical quantities in Sec. III.3, and the distribution functions of the total and staggered scalar chiralities in Sec. III.4. Summary and discussion of the results are given in Sec. IV, including the relationship between this study and previous studies for skyrmions in centrosymmetric metallic magnets. We show some of the details of the Ewald sum method in Appendix A, the data of the temperature dependences of the physical quantities in the single- $q$ phase in Appendix B, the definition and the temperature dependence of the $Z_{3}$ -symmetry breaking parameter in Appendix C, the histograms of the $Z_{3}$ -symmetry breaking parameter in each phase in Appendix D, the derivation of the four-body interaction induced by thermal fluctuations in Appendix E, and the evaluation of the small parameter $(J_{sd}/\epsilon_{\rm F})^{2}$ from the experimental data in Appendix F.

II Model and method

We consider the classical Heisenberg model on the 2D triangular lattice interacting with the long-range RKKY interaction. Our Hamiltonian is given by

	$\displaystyle H$	$\displaystyle=$	$\displaystyle-\sum_{i,j}J_{ij}\bm{S}_{i}\cdot\bm{S}_{j}-h\sum_{i=1}^{N}S_{i}^{z},$		(1)
	$\displaystyle J_{ij}$	$\displaystyle=$	$\displaystyle-J_{0}\left(\frac{\cos(2k_{\rm F}r_{ij})}{r_{ij}^{3}}-\frac{\sin(2k_{\rm F}r_{ij})}{2k_{\rm F}r_{ij}^{4}}\right),$		(2)

where $\bm{S}_{i}=(S_{i}^{x},S_{i}^{y},S_{i}^{z})$ is a three-component unit vector at the site $i$ , $J_{ij}$ is the RKKY interaction between the sites $i$ and $j$ with the energy scale $J_{0}>0$ ( $J_{0}$ is proportional to $\epsilon_{\rm F}(J/\epsilon_{\rm F})^{2}$ ), $r_{ij}=|\bm{r}_{i}-\bm{r}_{j}|$ is the distance between the sites $i$ and $j$ in units of the lattice constant $a$ , $k_{\rm F}$ is the Fermi wavenumber in units of $a^{-1}$ , and $h$ is the magnetic-field intensity. The sum is taken over all spin pairs on the 2D triangular lattice, and the total number of lattice sites is $N=L\times L$ . The dimensionless temperature and magnetic field are defined by $\tilde{T}=k_{\rm B}T/J_{0}$ and $\tilde{h}=h/J_{0}$ , and hereafter we denote $\tilde{T}$ and $\tilde{h}$ simply as $T$ and $h$ . Note that the RKKY interaction of Eq. (2) assumes the spherical Fermi surface of the 3D electron gas.

Refer to caption — Figure 1: The ground-state spin structures of the classical Heisenberg model on the 2D triangular lattice of (a) the RKKY model versus $k_{\rm F}$ , and of (c) the $J_{1}$ - $J_{3}$ model with the ferromagnetic nearest-neighbor coupling $J_{1}>0$ and the antiferromagnetic third-neighbor coupling $J_{3}<0$ versus $J_{1}/|J_{3}|$ . “ $120^{\circ}$ ” means the commensurate $120^{\circ}$ spin state with the three-sublattice periodicity, while “incommensurate” means the incommensurate spiral states. (b) The Fermi wavenumber $k_{\rm F}$ dependence of the ordering wavevector $|\bm{q}^{*}|$ around the incommensurate region of $2\pi/2.98<k_{\rm F}\leq 2\pi/2.58$ , which is obtained from the Ewald RKKY potential for $L=32$ .

Since the spin $\bm{S}_{i}$ is classical here, the ground state of the model in zero field can be obtained by the Fourier transform of $J_{ij}$ , which yields the spiral order for certain ranges of $k_{\rm F}$ with the incommensurate ordering wavevectors. Fig. 1 (a) shows the ground-state spin structures of the model as compared with those of the short-range $J_{1}$ - $J_{3}$ model of Ref. Okubo et al. (2012). In addition to the ferromagnetic and the commensurate $120^{\circ}$ -type orders, the incommensurate spiral order is realized for finite ranges of $k_{\rm F}$ , similarly to the short-range $J_{1}$ - $J_{3}$ case. We show in Fig. 1 (b) the Fermi wavenumber $k_{\rm F}$ dependence of the ordering wavenumber of the incommensurate spiral, $|\bm{q}^{*}|$ , which is related to the sizes of skyrmions, i.e., the skyrmion diameter being given by $\frac{4\pi}{\sqrt{3}|\bm{q}^{*}|}$ . As can be seen from Fig. 1 (b), the ordering wavenumbers cover a wide range of $|\bm{q}|$ spanning continuously from the short-wavelength $120^{\circ}$ structure to the long-wavelength uniform state. In this study, we concentrate on the incommensurate case and mostly set $k_{\rm F}=2\pi/2.71$ . The ground-state wavenumbers are three-fold degenerate as $q_{1}^{*},q_{2}^{*}$ and $q_{3}^{*}$ , reflecting the three-fold rotational symmetry of the triangular lattice.

Here, we comment on the relation of the present RKKY interaction (2) with the spin-quadratic interaction employed in the $q$ -space Hamiltonian of Refs. Hayami et al. (2017); Hayami and Motome (2019); Hayami and Yambe (2020); Hayami and Motome (2021a, b); Yambe and Hayami (2021), i.e., $\sum_{\bm{q}=\bm{q}^{*}}J_{\bm{q}}|\bm{S}_{\bm{q}}|^{2}$ , where the wavevectors $q_{1}^{*}\sim q_{3}^{*}$ are considered to be determined by the Fermi-surface effect (nesting). When such a quadratic interaction is rewritten in real space, one gets $\sum_{i<j}\sum_{\bm{q}=\bm{q}^{*}}\cos(\bm{q}\cdot(\bm{r}_{i}-\bm{r}_{j}))\bm{S}_{i}\cdot\bm{S}_{j}$ . This quadratic interaction oscillates in sign depending on the distance, just as in the RKKY interaction (2), yielding frustration in the standard sense, like the frustration in metallic spin glasses Mézard et al. (1987); Mydosh (1993); Kawamura and Taniguchi (2015). By contrast, it does not decay spatially even in the long-distance limit and corresponds to the mean-field approximation. Such an artificial feature arises due to the negligence of wavevector fluctuations around $q^{*}$ . In the present study, we employ the “textbook” RKKY interaction, i.e., we take account of the spatial decay of the interaction caused by the wavevector fluctuations without assuming a specific Fermi-surface effect.

The lattice sizes studied are $L=32,48,64$ with periodic boundary conditions. To fully take account of the long-range nature of the RKKY interaction beyond the finite system size $L$ , we employ the Ewald sum method, which is a general method to treat the long-range interaction in numerical simulations Ewald (1921); Hansen (1973); Fuchizaki (1994); Ikeda and Kawamura (2008); Mitsumoto and Kawamura (2021). Details of the derivation of the Ewald potential for the RKKY interaction on the triangular lattice are given in Appendix A.

To equilibrate the system, we combine the Metropolis and over-relaxation methods. We employ the replica-exchange method at relatively high temperatures and the simulated annealing method at lower temperatures. We also employ the mixed-phase method Creutz et al. (1979) to determine a part of the phase boundary. Our unit Monte Carlo step (MCS) consists of one Metropolis sweep and $L$ times over-relaxation sweeps. We take $1.0\times 10^{5}$ MCS each for equilibration and computing thermal averages in the replica-exchange simulations. Errors of physical quantities are estimated over three independent runs.

III Monte Carlo results

III.1 Phase diagram

In this section, we present our simulation results. We first show in Fig. 2 the temperature $T$ versus the magnetic-field $h$ phase diagram of the model with $k_{\rm F}=2\pi/2.71$ . The triple- $q$ SkX phase is stabilized at finite temperatures and intermediate fields, in addition to the single- $q$ , double- $q$ , and $Z$ phases. The obtained phase diagram is similar to the one of the short-range $J_{1}$ - $J_{3}$ model Okubo et al. (2012), including the properties of each phase as will be shown in the following subsections and sections.

III.2 Spin configurations in the wavevector- and real-spaces

Next, we wish to characterize each phase appearing in the phase diagram of Fig.2 by presenting our MC results of the spin structure factors, and the spin and chirality configurations in real space. The spin structure factor for the spin $xy$ -components $S_{\perp}(\bm{q})$ and the one for the spin $z$ -component $S_{\parallel}(\bm{q})$ are defined by

	$\displaystyle S_{\perp}(\bm{q})$	$\displaystyle=$	$\displaystyle\frac{1}{N}\left\langle\left(\sum_{\mu=x,y}\left\|\sum_{i=1}^{N}S_{i}^{\mu}e^{-i\bm{q}\cdot\bm{r}_{i}}\right\|^{2}\right)^{1/2}\right\rangle,$		(3)
	$\displaystyle S_{\parallel}(\bm{q})$	$\displaystyle=$	$\displaystyle\frac{1}{N}\left\langle\left\|\sum_{i=1}^{N}S_{i}^{z}e^{-i\bm{q}\cdot\bm{r}_{i}}\right\|\right\rangle,$		(4)

where $\langle\cdots\rangle$ represents the thermal average.

The spin structure factors $S_{\perp}(\bm{q})$ and $S_{\parallel}(\bm{q})$ computed in the triple- $q$ SkX phase are shown in Figs. 3 (a) and (b), respectively. In the triple- $q$ SkX phase, both the spin $xy$ -components and $z$ -component are characterized by all three equivalent pairs of ordering wavevectors, $\pm\bm{q}_{1}^{*},~{}\pm\bm{q}_{2}^{*}$ and $\pm\bm{q}_{3}^{*}$ , keeping the $Z_{3}$ lattice-rotation symemtry. Note that, due to the $U(1)$ spin-rotation symmetry associated with the spin $xy$ -components, sharp-looking peaks in $S_{\perp}(\bm{q})$ might not be true Bragg peaks but rather quasi-Bragg peaks associated with algebraic spin correlations Mermin and Wagner (1966); Bruno (2001).

Spin and chirality configurations in real space in the triple- $q$ SkX phase are shown in Figs.3 (c)-(f). The local scalar chirality is defined for the upward (downward) triangles by $\chi_{\bigtriangleup(\bigtriangledown)}=\bm{S}_{i}\cdot\bm{S}_{j}\times\bm{S}_{k}~{}(i,j,k\in\bigtriangleup(\bigtriangledown))$ . Due to the underlying $Z_{2}$ spin-reflection symmetry, the triple- $q$ SkX and anti-SkX states are energetically degenerate, either of which is realized as a result of the spontaneous $Z_{2}$ -symemtry breaking, as shown in the spin configurations of Figs. 3 (c,e) and the chirality configurations of Figs. 3 (d,f). The SkX and anti-SkX states have net total chiralities of mutually opposite signs, giving rise to the topological Hall effect of mutually opposite signs.

The spin structure factors of the $Z$ phase are shown in Figs. 4 (a) and (b). While they look more or less similar to those of the triple- $q$ SkX phase, the intensities in $S_{\perp}(\bm{q})$ are rather broad with finite widths, not being even the quasi-Bragg peaks, in sharp contrast to the intensities in $S_{\parallel}(\bm{q})$ . Such sharp versus broad features in $S_{\parallel}({\bm{q}})$ and $S_{\perp}({\bm{q}})$ intensities might serve as an experimental indicator of the $Z$ phase in distinction with the triple- $q$ SkX phase where both $S_{\parallel}({\bm{q}})$ and $S_{\perp}({\bm{q}})$ exhibit sharp intensities.

As shown in the real-space spin and chirality configurations shown in Figs. 4 (c) and (d), the $Z$ phase is the random-domain state consisting of the SkX and anti-SkX domains. On lowering $T$ toward the transition point to the triple- $q$ SkX phase, the sizes of these random domains grow, eventually leading to the triple- $q$ SkX state characterized by the net total chirality. Such behaviors of the spin structure factors and the real-space spin and chirality configurations are quite similar to those found in the short-range $J_{1}$ - $J_{3}$ model Okubo et al. (2012). Note that in both cases of the short-range and long-range models the finite-size SkX and anti-SkX domains in the $Z$ phase would be subject to thermal fluctuations, while more macroscopic SkX and anti-SkX domains in the triple- $q$ phase would be more or less frozen. Such a difference might be distinguishable by appropriate experimental means, e.g., by Lorentz Transmission Electron Microscopy.

The spin structure factors and the real-space spin and chirality configurations in the single- $q$ phase are shown in Figs. 5 (a,b) and 5 (c,d), respectively. $S_{\perp}(\bm{q})$ is characterized by one pair of ordering wavevectors selected from three equivalent wavevector pairs as shown in Fig. 5 (a), spontaneously breaking the lattice-rotation symmetry, which corresponds to the spiral direction in real space as shown in Fig 5 (c). By contrast, $S_{\parallel}(\bm{q})$ exhibits only weak broad intensities at $\pm\bm{q}_{1}$ , $\pm\bm{q}_{2}$ and $\pm\bm{q}_{3}$ . As shown in Fig.5 (d), the scalar chirality in the single- $q$ phase exhibits a staggered order, i.e., the upward and downward triangles have mutually opposite signs of the local scalar chirality.

The spin structure factors and the real-space spin and chirality configurations in the double- $q$ phase are shown in Figs. 6 (a,b) and 6 (c,d), respectively. $S_{\perp}(\bm{q})$ in the double- $q$ phase is characterized by two pairs of ordering wavevectors. The state consists of the superposition of the two-spirals formed by the spin $xy$ -components running along different directions as shown in the real-space spin configuration of Fig. 6 (c). By contrast, the spin $z$ -component forms the linear spin-density wave running along the direction complementary to those of the $xy$ -components, and the corresponding one-pair peak is observed in $S_{\parallel}(\bm{q})$ as shown in Fig. 6 (b). The scalar chirality shown in Fig. 6 (d) also exhibits a linear density wave in the same direction as that of the spin $z$ -component.

III.3 Temperature dependences of several physical quantities

To further clarify the properties of each phase, we next investigate the temperature dependences of several physical quantities, including the specific heat $C$ , the uniform magnetic susceptibility $\chi_{m}$ , the total scalar chirality $\chi_{\rm tot}$ and the staggered scalar chirality $\chi_{\rm stg}$ . The specific heat and the susceptibility are defined by the energy and the magnetization fluctuations, respectively. The total and staggerd scalar chiralities, which can be regarded as the order parameters of the triple- $q$ SkX and single- $q$ phases, respectively, are defined by

	$\displaystyle\chi_{\rm tot}$	$\displaystyle=\frac{1}{2N}\left(\sum_{\bigtriangleup}\chi_{\bigtriangleup}+\sum_{\bigtriangledown}\chi_{\bigtriangledown}\right),$		(5)
	$\displaystyle\chi_{\rm stg}$	$\displaystyle=\frac{1}{2N}\left(\sum_{\bigtriangleup}\chi_{\bigtriangleup}-\sum_{\bigtriangledown}\chi_{\bigtriangledown}\right),$		(6)

where the summation $\sum_{\bigtriangleup}$ ( $\sum_{\bigtriangledown}$ ) runs over all upward (downward) triangles on the triangular lattice. In the following, we show our MC data obtained by fully equilibrated replica-exchange simulations at two typical magnetic fields, i.e., the intermediate field of $h=0.10$ crossing in the phase diagram the paramagnetic, $Z$ , and triple- $q$ SkX phases, and the high field of $h=0.20$ crossing in the phase diagram the paramagnetic, single- $q$ , and double- $q$ phases.

Figs. 7 (a)-(d) shows the data at the intermediate magnetic field of $h=0.10$ in the temperature range of $T\geq 0.035$ . As can be seen from Fig. 7 (a), the specific heat exhibits two peaks at $T_{\rm c1}^{\rm(mid)}=0.0419$ and at $T_{\rm c2}^{\rm(mid)}=0.0400$ . As shown in Fig. 7 (c), the total scalar chirality $\chi_{\rm tot}$ is suppressed at $T_{\rm c2}^{\rm(mid)}<T<T_{\rm c1}^{\rm(mid)}$ but is enhanced at $T<T_{c2}^{\rm(mid)}$ . By contrast, the staggered scalar chirality $\chi_{\rm stg}$ shown in Fig. 7 (d) is vanishing in the investigated temperature range. Such behaviors are consistent with the identification of $T_{\rm c1}^{\rm(mid)}$ and $T_{\rm c2}^{\rm(mid)}$ as the paramagnetic/ $Z$ and the $Z$ /triple- $q$ transition temperatures, respectively. Note that, on further lowering $T$ to $T\leq 0.035$ , the system exhibits another phase transition into the single- $q$ phase, as shown in the phase diagram of Fig. 2.

Figs. 8 (a)-(d) shows the data at the high magnetic field of $h=0.20$ in the temperature range of $T\geq 0.025$ . Although the specific field shown in Fig. 8 (a) exhibits only a single peak at $T_{\rm c1}^{\rm(high)}=0.0312$ , the staggered scalar chirality $\chi_{\rm stg}$ shown in Fig. 8 (d) is enhanced in the intermediate temperature range $T_{\rm c2}^{\rm(high)}=0.0282<T<T_{\rm c1}^{\rm(high)}$ , but is suppressed below $T_{\rm c2}^{\rm(high)}$ . In addition, the behaviors of the $Z_{3}$ -symmetry-breaking parameters both for the spin $xy$ and $z$ components, shown in Appendix C for each typical magnetic field, suggest that the $Z_{3}$ lattice-rotation symmetry is broken only for the spin $xy$ components but not for the $z$ component at $T_{\rm c2}^{\rm(high)}<T<T_{\rm c1}^{\rm(high)}$ , while it is broken both for the spin $xy$ and $z$ components at $T<T_{\rm c2}^{\rm(high)}$ . The total scalar chirality $\chi_{\rm tot}$ shown in Fig. 8 (c) tends to vanish at any temperature. These observations are consistent with the identification of $T_{\rm c1}^{\rm(high)}$ and $T_{\rm c2}^{\rm(high)}$ as the paramagnetic/single- $q$ and the single- $q$ /double- $q$ transition temperatures, respectively. The magnetic susceptibility $\chi_{m}$ shown in Fig. 8 (b) exhibits a behavior more or less similar to the specific heat. Note that, on further lowering the temperature to $T\leq 0.025$ , the system exhibits another phase transition into the single- $q$ phase, as shown in the phase diagram of Fig. 2.

The corresponding data at the low field of $h=0.05$ , crossing in the phase diagram the paramagnetic and single- $q$ phases, are shown in Appendix B.

III.4 Distribution functions of the chirality

Very recently, an intriguing RSB phenomenon was observed in the triple- $q$ SkX phase and the double- $q$ phase in the 3D RKKY Heisenberg model on a stacked-triangular lattice where the multiple- $q$ states macroscopically coexist with the single- $q$ state Mitsumoto and Kawamura (2021). Under such circumstances, we wish to examine in this subsection whether a similar RSB occurs or not in the present 2D RKKY model. For this purpose, we study the distribution functions of the total and scalar chiralities in each phase.

In Fig. 9 (a), we show the distribution function of the total scalar chirality $P(\chi_{\rm tot})$ at $T=0.035$ and $h=0.1$ lying in the triple- $q$ SkX phase. Reflecting the $Z_{2}$ chiral degeneracy, there are positive and negative peaks of equal heights, each corresponding to the anti-SkX and SkX states, respectively, which grow and sharpen with $L$ . By contrast, no appreciable weight is found around $\chi_{\rm tot}=0$ corresponding to the single- $q$ state, in sharp contrast to the 3D case. Hence, no signature of the single- $q$ state coexisting with the triple- $q$ SkX (anti-SkX) state is observed, suggesting the absence of the RSB. Noting that the order parameter of the single- $q$ state is the staggered scalar chirality $\chi_{\rm stag}$ , we show in Fig. 9 (b) its distribution function $P(\chi_{\rm stag})$ . Only a sharp peak around zero corresponding to the triple- $q$ state is observed without any sign of the single- $q$ state, again indicating the absence of the RSB.

Likewise, there is no sign of the RSB in other phases, whose distribution functions of the total and staggered chiralities are shown in Figs. 9 (c)-(h). The staggered scalar chirality becomes finite only in the single- $q$ phase, and except in the single- $q$ phase itself, $P(\chi_{\rm stg})$ exhibits only a single peak around $\chi_{\rm stg}=0$ , indicating that the single- $q$ state coexists with neither the double- $q$ state nor the $Z$ state.

We also present the distribution functions of the $Z_{3}$ symmetry breaking order parameter in Appendix D. The distributions of the $Z_{3}$ symmetry breaking order parameter also indicate that the RSB does not occur in the 2D RKKY model.

IV Summary and Discussion

In summary, we have shown by extensive MC simulations on the 2D RKKY model without nesting, four-body interaction or magnetic anisotropy that the SkX (anti-SkX) phase is stabilized under magnetic fields only by the isotropic RKKY interaction. The SkX is stabilized by thermal fluctuations and frustration inherent to the RKKY interactions decaying as $1/r^{3}$ with the alternating signs. We have also shown the RSB does not occur in the 2D RKKY model, in sharp contrast to the 3D RKKY model.

We wish to discuss first why the RSB does not occur in the present 2D RKKY model, unlike the corresponding 3D RKKY model. In the fields of spin glasses and structural glasses where the RSB was first introduced and discussed, understanding of the RSB in finite dimensions has been regarded as a very hard and challenging problem and is still debated Binder, K. and Young, A. P. (1986); Mézard et al. (1987); Kawamura and Taniguchi (2015); Parisi, G. and Urbani, P. and Zamponi, F. (2020). In the case of glassy systems, while the occurrence of the RSB has been established in infinite dimensions, the situation is not necessarily clear in finite dimensions, though the general expectation might be that the RSB is more likely to arise in higher dimensions. Such a general tendency seems consistent with our present observation that the RSB is realized in the 3D RKKY system but not in the 2D RKKY system.

Here, we wish to give a further discussion on the issue of the existence/nonexistence of the RSB in the present RKKY model from the viewpoint of the underlying frustration effect, a crucial ingredient to realize the SkX state. In the 3D RKKY model studied in Ref.Mitsumoto and Kawamura (2021), not only the in-plane interaction but also the inter-plane interaction is frustrated due to the oscillating nature of the RKKY interaction. Although the ground state of the 3D RKKY model studied in Ref.Mitsumoto and Kawamura (2021) exhibits the ferromagnetic alignment along the inter-plane direction, apparently disguising the absence of inter-plane frustration, the ferromagnetic inter-plane spin alignment actually frustrates with the RKKY interaction which oscillates in sign, so that the 3D RKKY model is frustrated not only in the in-plane directions but also in the inter-plane direction. In that sense, the 3D RKKY model is more heavily frustrated than the 2D RKKY model, the latter possessing only the in-plane frustration.

The condition of the RSB is that the free energy difference between the macroscopically coexisting phases, which are the triple- $q$ SkX state and the single- $q$ state in the 3D RKKY model, stays only of $O(1)$ even in the thermodynamic limit. Such an unusual situation is expected to arise only in the presence of quite heavy frustration. Thus, the absence/presence of the RSB in the 2D/3D RKKY model might intuitively be understandable from the difference in the severity of the underlying frustration effect between the 2D/3D RKKY models.

In this connection, if one considers the 3D version of the short-range $J_{1}$ - $J_{3}$ model with the ferromagnetic inter-plane coupling, the SkX phase there is expected not to exhibit the RSB since the ferromagnetic inter-plane interaction does not bring about any additional frustration along the inter-plane direction. Indeed, the ordering properties of such a 3D model with the short-range interactions have been studied quite recently, indicating that the triple- $q$ SkX phase is realized, but not accompanying the RSB even in 3D Osamura, R. and Aoyama, K. and Mitsumoto, K. and Kawamura, H. . This observation highlights and confirms the importance of the severe frustration effect caused by the inter-plane RKKY interaction in realizing the RSB phenomenon.

Next, we discuss the possible implications of our results to the SkX formation mechanism. The nesting property of the Fermi surface, which has been supposed to be important in the SkX formation in previous studies Ozawa et al. (2017); Hayami et al. (2017); Hayami and Yambe (2020), might play a role in picking up particular wavevectors, but frustration alone can pick up appropriate ordering wavevectors even without nesting as our present model calculation has shown. In that sense, nesting plays only a secondary role in the SkX formation itself. As the system without nesting can exhibit a stable SkX phase, and the tuning of the interaction is usually necessary to realize the strong nesting, it would be important to check experimentally whether the candidate material really possesses strong nesting or not by independent measurements, e.g., the angle-resolved photoemission spectroscopy Inosov et al. (2009); Takahashi et al. (2020).

As discussed in Ref. Okubo et al. (2012) at the Ginzburg-Landau level, the effect of emergent many-body interactions generated by thermal effects is important in stabilizing the SkX state. In fact, as shown in Appendix F, such emergent many-body interactions can rigorously be derived in the form of the expansion in terms of the coarse-grained effective fields so long as $T$ is not too low compared with $T_{c}$ Fisher (1983); Kawamura (1988). We note that these many-body terms, including the biquadratic interaction with the positive coefficient, have the same forms as those discussed in Ref.Hayami et al. (2017) in the context of the electronic perturbation with respect to $J_{sd}/\epsilon_{\rm F}$ . Note that while the energy scale of the four-body (biquadratic) interaction of Ref.Hayami et al. (2017) is smaller than that of the quadratic RKKY interaction $J_{0}$ by the factor of $(J_{sd}/\epsilon_{\rm F})^{2}$ , the energy scale of the emergent four-body interaction arising from thermal fluctuations is of the same order as $J_{0}$ , not containing the small parameter $(J_{sd}/\epsilon_{\rm F})^{2}$ . With Gd₂PdSi₃ in mind, we estimate from the available experimental data the typical $(J_{sd}/\epsilon_{\rm F})^{2}$ value to be as small as $\sim 10^{-3}-10^{-4}$ : Details are given in Appendix F. Hence, at least at finite $T$ , the emergent many-body interaction originating from thermal fluctuations would play a dominant role in stabilizing the SkX.

Interesting questions are whether the SkX is stabilized at $T\rightarrow 0$ and its stabilization mechanism. It has been reported experimentally that the SkX state, once stabilized at finite temperatures, can be continued down to lower temperatures as a metastable state Oike, H. and Kikkawa, A. and Kanazawa, N. and Taguchi, Y. and Kawasaki, M. and Tokura, Y. and Kagawa, F. (2016); Karube, K. and White, J. S. and Reynolds, N. and Gavilano, J. L. and Oike, H. and Kikkawa, A. and Kagawa, F. and Tokunaga, Y. and Rønnow, H. M. and Tokura, Y. and Taguchi, Y. (2016). One possible mechanism to realize the stable $T=0$ SkX state might be the biquadratic interaction of electronic origin as discussed in Refs.Hayami et al. (2017). We note that quantum fluctuations of localized spins often play similar roles as thermal fluctuations in frustrated magnets, stabilizing the nontrivial multiple- $q$ state even at $T=0$ even when the classical ground state is of the single- $q$ type Marmorini and Momoi (2014). Again, such a term originating from quantum spin fluctuations does not contain the small parameter $(J_{sd}/\epsilon_{\rm F})^{2}$ and might play a role in the SkX stabilization at lower $T$ , especially when the spin quantum number $S$ is small. Although we limit our discussion so far to isotropic Heisenberg spins, real magnets possess some amount of magnetic anisotropy. Indeed, previous theoretical studies have indicated that the easy-axis-type anisotropy tends to stabilize the SkX state even at $T=0$ both for the classical Leonov and Mostovoy (2015); Hayami et al. (2016); Wang et al. (2020, 2021) and quantum Lohani, V. and Hickey, C. and Masell, J. and Rosch, A. (2019) systems. Thus, the magnetic anisotropy might actually play a major role in the SkX stabilization at $T=0$ .

Acknowledgements.

The authors would like to thank K. Aoyama for useful discussion. We are thankful to ISSP, the University of Tokyo, and YITP, Kyoto University, for providing us with CPU time. This work is supported by JSPS KAKENHI Grant No. JP17H06137.

Appendix A Ewald-sum method of the RKKY interaction on the two-dimensional triangular lattice

In this section, we explain some of the details of the application of the Ewald-sum method, widely used in computing the long-range interactions in periodic systems of finite sizes Ewald (1921); Hansen (1973); Fuchizaki (1994); Ikeda and Kawamura (2008); Mitsumoto and Kawamura (2021), to the RKKY Heisenberg model on the 2D triangular lattice. In fact, in the straightforward application of the Ewald method to the present 2D problem, we meet some computational difficulty due to poor numerical convergence. In order to circumvent this difficulty, we first consider the 3D system and then take the 2D limit, i.e., the limit of the ratio $c$ between the interlayer and intralayer distances tending to infinity at the end.

We begin with the Ewald potential for the RKKY interaction on the 3D stacked-triangular lattice employed in Ref.Mitsumoto and Kawamura (2021),

$\displaystyle J_{ij}^{\rm Ewald(3D)}$	$\displaystyle=$	$\displaystyle-J_{0}\sum_{\bm{\lambda}}J_{ij}(\bm{\lambda}),$	(7)
$\displaystyle J_{ij}(\bm{\lambda})$	$\displaystyle=$	$\displaystyle\frac{\cos(2k_{\rm F}\|\bm{r}_{ij}+L\bm{\lambda}\|)}{\|\bm{r}_{ij}+L\bm{\lambda}\|^{3}}$	(8)
	$\displaystyle-$	$\displaystyle\frac{\sin(2k_{\rm F}\|\bm{r}_{ij}+L\bm{\lambda}\|)}{2k_{\rm F}\|\bm{r}_{ij}+L\bm{\lambda}\|^{4}},$	(8)

where $L$ is linear size of the system, and $\bm{\lambda}=n_{a}\bm{a}+n_{b}\bm{b}+n_{c}\bm{c}$ with $\bm{a}=(1,0,0),~{}\bm{b}=(1/2,\sqrt{3}/2,0),~{}\bm{c}=(0,0,c)$ . The sum over $\bm{\lambda}$ in Eq. (7) runs over integers $n_{\mu}=-\infty,...,0,...\infty$ ( $\mu=a,b,c)$ , $L\bm{\lambda}$ mapping the original cell of the size $L^{3}$ to the image cell with exactly the same spin configuration as that in the original cell.

We divide the potential into two parts, a short-range contribution and a long-range contribution, by inserting the identity,

1=\frac{1}{\Gamma(\frac{3}{2})}\left[\Gamma\left(\frac{3}{2},\pi\frac{|\bm{r}_{ij}+L\bm{\lambda}|^{2}}{L^{2}}\right)+\gamma\left(\frac{3}{2},\pi\frac{|\bm{r}_{ij}+L\bm{\lambda}|^{2}}{L^{2}}\right)\right],

(9)

where $\Gamma(\alpha)=\int_{0}^{\infty}dte^{-t}t^{\alpha-1}$ is the gamma function, $\Gamma(\alpha,x)=\int_{x}^{\infty}dte^{-t}t^{\alpha-1}$ and $\gamma(\alpha,x)=\int_{0}^{x}dte^{-t}t^{\alpha-1}$ are the incomplete gamma functions. The short-range contribution is well-converged in the real space. To also make the long-range contribution well-converged, we transform it into the Fourier space. Then, Eq. (7) can be rewritten as Mitsumoto and Kawamura (2021),

$\displaystyle J_{ij}^{\rm Ewald(3D)}$	$\displaystyle=$	$\displaystyle-\frac{J_{0}}{\Gamma(\frac{3}{2})}\left(J_{ij}^{\rm short(3D)}+J_{ij}^{\rm long(3D)}\right)$	(10)
$\displaystyle J_{ij}^{\rm short(3D)}$	$\displaystyle=$	$\displaystyle\sum_{\bm{\lambda}}\Gamma\left(\frac{3}{2},\pi\frac{\|\bm{r}_{ij}+L\bm{\lambda}\|^{2}}{L^{2}}\right)\left[\frac{\cos(2k_{\rm F}\|\bm{r}_{ij}+L\bm{\lambda}\|)}{\|\bm{r}_{ij}+L\bm{\lambda}\|^{3}}-\frac{\sin(2k_{\rm F}\|\bm{r}_{ij}+L\bm{\lambda}\|)}{2k_{\rm F}\|\bm{r}_{ij}+L\bm{\lambda}\|^{4}}\right],$	(11)
$\displaystyle J_{ij}^{\rm long(3D)}$	$\displaystyle=$	$\displaystyle\frac{\pi^{\frac{3}{2}}}{\sqrt{3}cL^{3}}\sum_{\bm{h}}\frac{e^{2\pi i\frac{\bm{h}}{L}}\cdot\bm{r}_{ij}}{h}\Bigg{\{}h_{+}E_{1}(\pi h_{+}^{2})-h_{-}E_{1}(\pi h_{-}^{2})$	(12)
	$\displaystyle-$	$\displaystyle\frac{1}{2k_{\rm F}L}\left[\pi h_{+}^{2}E_{1}(\pi h_{+}^{2})-\pi h_{-}^{2}E_{1}(\pi h_{-}^{2})-(e^{-\pi h_{+}^{2}}-e^{-\pi h_{-}^{2}})\right]\Bigg{\}},$	(12)

where $\bm{h}=m_{a}\bm{k}_{a}+m_{b}\bm{k}_{b}+m_{c}\bm{k}_{c}$ with $\bm{k}_{a}=(1,-1/\sqrt{3},0),~{}\bm{k}_{b}=(0,2/\sqrt{3},0),~{}\bm{k}_{c}=(0,0,1/c)$ , $h=|\bm{h}|$ , $h_{\pm}=h\pm Lk_{\rm F}/\pi$ , and $E_{1}(x)=\int_{x}^{\infty}dte^{-t}/t$ is the exponential integral function. The sum $\sum_{\bm{h}}$ runs over all integers $m_{\mu}=-\infty,...,0,...\infty$ ( $\mu=a,b,c$ ).

Then, we take the 2D limit, i.e., $c\to\infty$ . The Ewald potential can be written as

J_{ij}^{\rm Ewald(2D)}=-\frac{J_{0}}{\Gamma(\frac{3}{2})}\left(J_{ij}^{\rm short(2D)}+J_{ij}^{\rm long(2D)}\right).

(13)

The short-range contribution in 2D has the same form as that in 3D,

J_{ij}^{\rm short(2D)}=J_{ij}^{\rm short(3D)},

(14)

while $\bm{\lambda}$ is restricted to the $a$ - $b$ plane. In case of the long-range contribution, the summation in the $c$ -direction is replaced by the integral,

\frac{1}{c}\sum_{\bm{h}}=\frac{1}{c}\sum_{\bm{h}_{ab}}\sum_{\bm{h}_{c}}\underset{c\to\infty}{\to}\sum_{\bm{h}_{ab}}\int_{-\infty}^{\infty}dh_{c},

(15)

where $\bm{h}_{ab}=m_{a}\bm{k}_{a}+m_{b}\bm{k}_{b}$ and $\bm{h}_{c}=m_{c}\bm{k}_{c}$ . The long-range term $J_{ij}^{\rm long}$ then becomes

$\displaystyle J_{ij}^{\rm long(2D)}$	$\displaystyle=\frac{\pi^{\frac{3}{2}}}{\sqrt{3}L^{3}}\sum_{\bm{h}_{ab}}e^{2\pi i\frac{\bm{h}}{L}\cdot\bm{r}_{ij}}\int_{-\infty}^{\infty}dh_{c}$	(16)
	$\displaystyle\times\frac{1}{h}\Bigg{\{}h_{+}E_{1}(\pi h_{+}^{2})-h_{-}E_{1}(\pi h_{-}^{2})$
	$\displaystyle-\frac{1}{2k_{\rm F}L}\Big{[}\pi h_{+}^{2}E_{1}(\pi h_{+}^{2})-\pi h_{-}^{2}E_{1}(\pi h_{-}^{2})$
	$\displaystyle-(e^{-\pi h_{+}^{2}}-e^{-\pi h_{-}^{2}})\Big{]}\Bigg{\}}.$

The numerical evaluation of the right-hand side of Eq. (16) requires the evaluation of the double integral ( $E_{1}$ in the integrand contains an integral). In fact, the $\int_{-\infty}^{\infty}dh_{c}(h_{-}/h)E_{1}(\pi h_{-}^{2})$ term has an $x\log x$ -type singularity around $h_{-}=0$ . To properly take account of this singularity in the numerical integration for $h_{c}$ by using Simpson’s rule, we take a finer mesh of $10^{-5}$ around $h_{-}=0$ , while we take $10^{-4}$ in other intervals.

Appendix B Monte Carlo results of the temperature dependences of several physical quantities in single- $q$ phase

In Figs. 10 (a)-(d), we show the data at the low magnetic field of $h=0.05$ in the temperature range of $T\geq 0.037$ . The specific heat $C$ shown in Fig. 10 (a) exhibits a sharp peak at $T_{\rm c}^{\rm(low)}=0.0414$ , corresponding to the phase transition between the paramagnetic and single- $q$ phases. The magnetic susceptibility $\chi_{m}$ shown in Fig. 10 (b) exhibits a behavior more or less similar to that of the specific heat. As can be seen from Fig. 10 (d), the staggered scalar chirality $\chi_{\rm stg}$ is enhanced in the single- $q$ phase, while the total scalar chirality $\chi_{\rm tot}$ is not as shown in Fig. 10 (c).

Appendix C The $Z_{3}$ -symmetry-breaking parameters

The $Z_{3}$ -symmetry-breaking parameters for the spin $xy$ -components and the $z$ -component are defined by

	$\displaystyle\bm{\rho}_{3\perp}$	$\displaystyle=s_{\perp}(\bm{q}_{1}^{})\bm{e}_{1}+s_{\perp}(\bm{q}_{2}^{})\bm{e}_{2}+s_{\perp}(\bm{q}_{3}^{*})\bm{e}_{3},$			(17)
	$\displaystyle\bm{\rho}_{3\parallel}$	$\displaystyle=s_{\parallel}(\bm{q}_{1}^{})\bm{e}_{1}+s_{\parallel}(\bm{q}_{2}^{})\bm{e}_{2}+s_{\parallel}(\bm{q}_{3}^{*})\bm{e}_{3},$			(18)

where $\bm{e}_{1}=(1,0)$ , $\bm{e}_{2}=(-1/2,\sqrt{3}/2)$ and $\bm{e}_{3}=(-1/2,-\sqrt{3}/2)$ are the unit vectors, $s_{\perp}(\bm{q})$ and $s_{\parallel}(\bm{q})$ are the instantaneous spin structure factors given by Eqs. (3) and (4) without the thermal average. These order parameters $\bm{\rho}_{3\perp}$ and $\bm{\rho}_{3\parallel}$ are two-component vectors defined in the two-dimensional triangular order-parameter space and represent the extent of the $Z_{3}$ -symmetry breaking. Hereafter, we call them the perpendicular and parallel $Z_{3}$ -breaking parameters. $\bm{\rho}_{3\parallel}$ is regarded as the order parameter of the double- $q$ phase because the $Z_{3}$ symmetry is not broken for the spin $z$ -component in the other phases.

Figures 11 (a)-(f) show the $Z_{3}$ -breaking parameters at three typical magnetic fields: i.e., the intermediate field of $h=0.10$ crossing in the phase diagram the paramagnetic, $Z$ , and triple- $q$ SkX phases, the high field of $h=0.20$ crossing in the phase diagram the paramagnetic, single- $q$ , and double- $q$ phases, and the low field of $h=0.05$ crossing in the phase diagram the paramagnetic and single- $q$ phases. The perpendicular $Z_{3}$ -breaking parameter $|\bm{\rho}_{3\perp}|$ is enhanced in both the single- $q$ and double- $q$ phases, while the parallel $Z_{3}$ -symmetry breaking parameter $|\bm{\rho}_{3\parallel}|$ is enhanced only in the double- $q$ phase.

Appendix D Distribution functions of the $Z_{3}$ -breaking parameter

In this section, we show our MC results of the distributions of the $Z_{3}$ -breaking parameter in various phases. Attention is paid to the issue of the existence/nonexistence of the RSB recently identified in the 3D RKKY system where the single- $q$ state macroscopically coexists with the triple- $q$ SkX state or the double- $q$ state Mitsumoto and Kawamura (2021).

The non-RSB feature can also be seen in the histograms of the perpendicular and parallel $Z_{3}$ -breaking parameters, $P(\bm{\rho}_{3\perp})$ and $P(\bm{\rho}_{3\parallel})$ , shown in Fig. 12. The single- $q$ state is characterized by the intensities in $P(\bm{\rho}_{3\perp})$ located at the corners of the triangular region, as shown in the left panel of Fig. 12 (a). As can be seen from Figs. 12 (b)-(d), no such intensity is detected in either the double- $q$ , triple- $q$ , or $Z$ phases.

Appendix E Derivation of the effective many-body interactions due to thermal fluctuations

In this section, we rigorously derive from the microscopic Heisenberg spin Hamiltonian the effective many-body interactions due to thermal fluctuations by performing a kind of coarse graining procedure Fisher (1983); Kawamura (1988). We consider the $n$ -component classical vector spins $\bm{S}_{i}=(S_{i}^{1},S_{i}^{2},...,S_{i}^{n})~{}(i=1,2,...,N)$ normalized as $|\bm{S}_{i}|=1$ on an arbitrary lattice in general $d$ dimensions, interacting via only the two-body interactions $J_{ij}$ whose energy scale is $J_{0}$ . It includes the RKKY isotropic Heisenberg model treated in the main text as a special case. The Hamiltonian under the site-dependent generalized magnetic field $\bm{h}_{i}$ is given by

\beta H=-\frac{1}{2}\sum_{i=1}^{N}\sum_{j\neq i}J_{ij}\bm{S}_{i}\cdot\bm{S}_{j}-\sum_{i=1}^{N}\bm{h}_{i}\cdot\bm{S}_{i}.

(19)

Introducing the dimensionless matrix $\tilde{J}_{ij}=P_{0}\delta_{ij}+J_{ij}/J_{0}$ , where $P_{0}>0$ is a positive constant sufficiently large to make $\tilde{J}_{ij}$ positive definite, one can apply the Hubbard-Stratonovich transformation,

$\displaystyle Z$	$\displaystyle=$	$\displaystyle\underset{\{\bm{S_{i}}\}}{\rm Tr}\exp\left(\frac{\beta J_{0}}{2}\sum_{i=1}^{N}\sum_{j=1}^{N}\tilde{J}_{ij}\bm{S}_{i}\cdot\bm{S}_{j}+\beta\sum_{i=1}^{N}\bm{h}_{i}\bm{S}_{i}\right)$	(20)
	$\displaystyle=$	$\displaystyle\frac{1}{A}\int\prod_{i=1}^{N}(d\bm{\phi}_{i})\underset{\{\bm{S_{i}}\}}{\rm Tr}\exp\left(-\frac{\beta J_{0}}{2}\sum_{i=1}^{N}\sum_{j=1}^{N}(\tilde{J}^{-1})_{ij}\bm{\phi}_{i}\cdot\bm{\phi}_{j}+\sum_{i=1}^{N}(\beta J_{0}\bm{\phi}_{i}+\beta\bm{h}_{i})\cdot\bm{S}_{i}\right)$
	$\displaystyle=$	$\displaystyle\frac{1}{A}\int\prod_{i=1}^{N}(d\bm{\phi}_{i})e^{-\beta H_{\rm eff}},$

where $\beta=1/T$ is the inverse temperature, and $H_{\rm eff}$ is the effective Hamiltonian given by,

	$\displaystyle H_{\rm eff}$	$\displaystyle=J_{0}\Bigg{[}\frac{1}{2}\sum_{i=1}^{N}\sum_{j=1}^{N}(\tilde{J}^{-1})_{ij}\bm{\phi}_{i}\cdot\bm{\phi}_{j}$
		$\displaystyle-\frac{1}{\beta J_{0}}\log\underset{\{\bm{S_{i}}\}}{\rm Tr}\prod_{i=1}^{N}e^{(\beta J_{0}\bm{\phi}_{i}+\beta\bm{h}_{i})\cdot\bm{S}_{i}}\Bigg{]},$		(21)

with $A=((2\pi/\beta J_{0})^{N}/|{\rm det}\tilde{J}^{-1}|)^{n/2}$ .

The auxiliary field $\bm{\phi}_{i}$ and the $n$ -component spin $\bm{S}_{i}$ are related as,

\langle S_{i}^{\mu}\rangle=\frac{\partial\log Z}{\partial(\beta h_{i}^{\mu})}=\sum_{j=1}^{N}(\tilde{J}^{-1})_{ij}\langle\phi_{j}^{\mu}\rangle,

(22)

which can easily be confirmed by integrating by parts after shifting the derivative with respect to $\beta h_{i}^{\nu}$ to the one with respect to $\beta J_{0}\phi_{i}^{\nu}$ . Since the field $\bm{\phi}_{i}$ is of unconstrained length, Eq. (22) suggests that $\bm{\phi}_{i}$ can be regarded as the coarse-grained spin variable.

We consider below the zero field case for simplicity. The trace in the second term of Eq. (21) can be taken as

\log\underset{\{\bm{S_{i}}\}}{\rm Tr}\prod_{i=1}^{N}e^{\beta J_{0}\bm{\phi}_{i}\cdot\bm{S}_{i}}=\sum_{i=1}^{N}\log\left(\frac{(2\pi)^{\frac{n}{2}}I_{\frac{n}{2}-1}(\beta J_{0}\phi_{i})}{(\beta J_{0}\phi_{i})^{\frac{n}{2}-1}}\right),

(23)

where

I_{\nu}(x)=\sum_{m=0}^{\infty}\frac{1}{m!\Gamma(m+\nu+1)}\left(\frac{x}{2}\right)^{2m+\nu},

(24)

is the modified Bessel function and $\phi_{i}$ is the absolute value of $\bm{\phi}_{i}$ . We then expand the effective Hamiltonian up to the fourth order in $\bm{\phi}_{i}$ to get,

	$\displaystyle H_{\rm eff}$	$\displaystyle=J_{0}\Bigg{[}\frac{1}{2}\sum_{i=1}^{N}\sum_{j=1}^{N}(\tilde{J}^{-1})_{ij}\bm{\phi}_{i}\cdot\bm{\phi}_{j}$
		$\displaystyle-\frac{\beta J_{0}}{2n}\sum_{i=1}^{N}\phi_{i}^{2}+\frac{(\beta J_{0})^{3}}{4n^{2}(n+2)}\sum_{i=1}^{N}\phi_{i}^{4}\Bigg{]}.$		(25)

Using the translational symmetry, the effective Hamiltonian can be rewritten as

	$\displaystyle H_{\rm eff}$	$\displaystyle=J_{0}\Bigg{[}\frac{1}{2}\sum_{\bm{q}}\left(\tilde{J}^{-1}_{\bm{q}}-\frac{\beta J_{0}}{n}\right)\|\bm{\phi}_{\bm{q}}\|^{2}$
		$\displaystyle+\frac{(\beta J_{0})^{3}}{4n^{2}(n+2)}\sum_{\{\bm{q}\}}{}^{\prime}(\bm{\phi}_{\bm{q}_{1}}\cdot\bm{\phi}_{\bm{q}_{2}})(\bm{\phi}_{\bm{q}_{3}}\cdot\bm{\phi}_{\bm{q}_{4}})\Bigg{]},$		(26)

where $\bm{\phi}_{\bm{q}}$ and $\tilde{J}^{-1}_{\bm{q}}$ the Fourier transforms of $\bm{\phi}_{i}$ and $(\tilde{J}^{-1})_{ij}$ ,

	$\displaystyle\bm{\phi}_{\bm{q}}$	$\displaystyle=$	$\displaystyle\frac{1}{N}\sum_{i=1}^{N}\bm{\phi}_{i}e^{-i\bm{q}\cdot\bm{r}_{i}},$		(27)
	$\displaystyle\tilde{J}^{-1}_{\bm{q}}$	$\displaystyle=$	$\displaystyle\sum_{j=1}^{N}(\tilde{J}^{-1})_{ij}e^{-i\bm{q}\cdot\Delta\bm{r}_{ij}},$		(28)

and $\sum_{\{\bm{q}\}}{}^{\prime}$ represents the restricted summation under the constraint $\bm{q}_{1}+\bm{q}_{2}+\bm{q}_{3}+\bm{q}_{4}=\bm{0}$ . Note that the prefactor of the fourth-order term is positive. Indeed, the summation includes the $q$ -space biquadratic interaction term with the positive coefficient which is of exactly the same form as the one used in Ref. Hayami et al. (2017) if the relevant $\bm{q}$ ’s are restricted to the ordering wavevectors $\bm{q}^{*}$ . Near the critical point $T_{c}$ , the temperature is of order $J_{0}$ , and $\beta J_{0}$ is of $O(1)$ . Therefore, the many-body interactions due to thermal fluctuations have the same energy scale as the original two-body interaction, not containing the small parameter such as $(J_{sd}/\epsilon_{\rm F})$ . Meanwhile, as the temperature is further lowered away from $T_{c}$ , $\beta J_{0}$ becomes large, and the expansion of the effective Hamiltonian eventually breaks down.

Appendix F Evaluation of the small parameter $(J_{sd}/\epsilon_{\rm F})^{2}$ from the experimental data

In this section, we numerically evaluate the typical value of the small parameter $(J_{sd}/\epsilon_{\rm F})^{2}$ from the available experimental data on Gd₂PdSi₃, a centrosymmetric metallic magnet exhibiting the SkX state in applied magnetic fields. For this purpose, we first evaluate the Fermi energy $\epsilon_{F}$ from the angle-resolved photoemission spectroscopy data and the energy scale of the RKKY interaction $J_{0}$ from the high-temperature susceptibility data. Then, $J_{sd}$ is related to $J_{0}$ and $\epsilon_{F}$ as Yosida (1957)

J_{0}\approx 9\pi\epsilon_{\rm F}\left(\frac{J_{sd}}{\epsilon_{\rm F}}\right)^{2},

(29)

where the half-filling condition is assumed for simplicity.

The Fermi energy $\epsilon_{F}$ of Gd₂PdSi₃ might be estimated from the angle-resolved photoemission spectroscopy measurements. The energy of the Fermi level measured from the bottom of the conduction band was reported to be $\epsilon_{\rm F}=0.50$ eV $=5.8\times 10^{3}$ K Inosov et al. (2009).

We estimate $J_{0}$ of Gd₂PdSi₃ by comparing the Curie-Weiss temperature determined experimentally from the high-temperature susceptibility measurements and the one determined theoretically from the high-temperature expansion of the 2D RKKY model. At high temperatures, the inverse magnetic susceptibility exhibits a linear dependence on the temperature $T$ as,

\chi^{-1}\propto(T-\Theta_{\rm p}),

(30)

where $\Theta_{p}$ is the Curie-Weiss temperature. According to Ref. Mallik et al. (1998), $\Theta_{p}$ of Gd₂PdSi₃ is positive, i.e., ferromagnetic, and its value was estimated to be $\sim$ 19K. Approximating Gd₂PdSi₃ by the 2D RKKY model on the triangular lattice described by Eq.(1) of the main text, we derive the high-temperature form of the inverse susceptibility by means of the standard high-temperature expansion. Expansion to the first nontrivial order yields,

\chi^{-1}=\frac{3}{4}\left[T-\frac{1}{3}\sum_{j=1}^{N}J_{ij}\right].

(31)

From Eqs. (30) and (31), we get

\Theta_{\rm p}(J_{0},k_{F})=\frac{1}{3}\sum_{j=1}^{N}J_{ij}(J_{0},k_{F}).

(32)

Note that the RKKY interaction $J_{ij}$ , and consequently $\Theta_{p}$ , depend not only on $J_{0}$ but also on $k_{F}$ . Generally speaking, the sign of $\Theta_{p}$ could be either positive (ferromagnetic) or negative (antiferromagnetic) depending on the $k_{F}$ -value. In view of the experimental observation of the positive $\Theta_{p}$ for Gd₂PdSi₃, we consider here the $k_{F}$ -region corresponding to the incommensurate spiral order with the positive $\Theta_{p}$ in the ground state. Referring to the ground-state phase diagram shown in Fig.1 of the main text, we choose here two candidate $k_{F}$ values, i.e., $k_{\rm F}=2\pi/2.71$ and $2\pi/2.61$ , the former being exactly the $k_{F}$ -value employed in the MC simulation of the main text, and numerically estimate the right-hand side of Eq. (32) based on the Ewald-sum method with $L=32$ (the $L$ -dependence turns out to be almost negligible). From the relation $\Theta_{p}=19$ K, we numerically estimate $J_{0}$ to be 250K and 59K, and by using Eq. (29) the small parameter $(J_{sd}/\epsilon_{\rm F})^{2}$ to be $1.5\times 10^{-3}$ and $3.6\times 10^{-4}$ for $k_{\rm F}=2\pi/2.71$ and $2\pi/2.61$ , respectively. Hence, the parameter $(J_{sd}/\epsilon_{\rm F})^{2}$ is likely to be quite small in the metallic magnet Gd₂PdSi₃, around $10^{-3}-10^{-4}$ . This observation indicates that the fourth-body interaction arising from the higher-order perturbation of Ref. Hayami et al. (2017) would be smaller than the quadratic RKKY interaction by the factor of $10^{-3}-10^{-4}$ in typical weak-coupling metals.

References

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–919 (2009).
Neubauer et al. (2009) A. Neubauer, C. Pfleiderer, B. Binz, A. Rosch, R. Ritz, P. G. Niklowitz, and P. Böni, “Topological Hall effect in the ${A}$ phase of MnSi,” Phys. Rev. Lett. 102, 186602 (2009).
Münzer et al. (2010) W. Münzer, A. Neubauer, T. Adams, S. Mühlbauer, C. Franz, F. Jonietz, R. Georgii, P. Böni, B. Pedersen, M. Schmidt, A. Rosch, and C. Pfleiderer, “Skyrmion lattice in the doped semiconductor Fe_1-xCo_xSi,” Phys. Rev. B 81, 041203 (2010).
Yu et al. (2010) X. Z. Yu, Y. Onose, N. Kanazawa, J. H. Park, J. H. Han, Y. Matsui, N. Nagaosa, and Y. Tokura, “Real-space observation of a two-dimensional skyrmion crystal,” Nature (London) 465, 901–904 (2010).
Yu et al. (2011) X. Z. Yu, N. Kanazawa, Y. Onose, K. Kimoto, W. Z. Zhang, S. Ishiwata, Y. Matsui, and Y. Tokura, “Near room-temperature formation of a skyrmion crystal in thin-films of the helimagnet FeGe,” Nat. Mater. 10, 106–109 (2011).
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).
Kurumaji et al. (2019) T. Kurumaji, T. Nakajima, M. Hirschberger, A. Kikkawa, Y. Yamasaki, H. Sagayama, H. Nakao, Y. Taguchi, T.-h. Arima, and Y. Tokura, “Skyrmion lattice with a giant topological Hall effect in a frustrated triangular-lattice magnet,” Science 365, 914–918 (2019).
Takahashi et al. (2020) H. Takahashi, K. Aono, Y. Nambu, R. Kiyanagi, T. Nomoto, M. Sakano, K. Ishizaka, R. Arita, and S. Ishiwata, “Competing spin modulations in the magnetically frustrated semimetal EuCuSb,” Phys. Rev. B 102, 174425 (2020).
Leonov and Mostovoy (2015) A. O. Leonov and M. Mostovoy, “Multiply periodic states and isolated skyrmions in an anisotropic frustrated magnet,” Nat. commun. 6, 8275 (2015).
Hayami et al. (2016) S. Hayami, S.-Z. Lin, and C. D. Batista, “Bubble and skyrmion crystals in frustrated magnets with easy-axis anisotropy,” Phys. Rev. B 93, 184413 (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 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).
Lin and Batista (2018) S.-Z. Lin and Cristian D. Batista, “Face centered cubic and hexagonal close packed skyrmion crystals in centrosymmetric magnets,” Phys. Rev. Lett. 120, 077202 (2018).
Hayami and Motome (2019) S. Hayami and Y. Motome, “Effect of magnetic anisotropy on skyrmions with a high topological number in itinerant magnets,” Phys. Rev. B 99, 094420 (2019).
Wang et al. (2020) Z. Wang, Y. Su, S.-Z. Lin, and C. D. Batista, “Skyrmion crystal from RKKY interaction mediated by 2D electron gas,” Phys. Rev. Lett. 124, 207201 (2020).
Hayami and Yambe (2020) S. Hayami and R. Yambe, “Degeneracy lifting of néel, bloch, and anti-skyrmion crystals in centrosymmetric tetragonal systems,” J. Phys. Soc. Jpn. 89, 103702 (2020), https://doi.org/10.7566/JPSJ.89.103702 .
Hayami and Motome (2021a) S. Hayami and Y. Motome, “Noncoplanar multiple- $q$ spin textures by itinerant frustration: Effects of single-ion anisotropy and bond-dependent anisotropy,” Phys. Rev. B 103, 054422 (2021a).
Hayami and Motome (2021b) S. Hayami and Y. Motome, “Square skyrmion crystal in centrosymmetric itinerant magnets,” Phys. Rev. B 103, 024439 (2021b).
Yambe and Hayami (2021) R. Yambe and S. Hayami, “Skyrmion crystals in centrosymmetric itinerant magnets without horizontal mirror plane,” Sci. Rep. 11, 11184 (2021).
Wang et al. (2021) Z. Wang, Y. Su, S.-Z. Lin, and C. D. Batista, “Meron, skyrmion, and vortex crystals in centrosymmetric tetragonal magnets,” Phys. Rev. B 103, 104408 (2021).
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 Gd₂PdSi₃,” Phys. Rev. B 60, 12162–12165 (1999).
Ruderman and Kittel (1954) M. A. Ruderman and C. Kittel, “Indirect exchange coupling of nuclear magnetic moments by conduction electrons,” Phys. Rev. 96, 99–102 (1954).
Kasuya (1956) T. Kasuya, “A theory of metallic ferro-and antiferromagnetism on Zener’s model,” Prog. Theor. Phys. 16, 45–57 (1956).
Yosida (1957) K. Yosida, “Magnetic properties of Cu-Mn alloys,” Phys. Rev. 106, 893–898 (1957).
Mitsumoto and Kawamura (2021) K. Mitsumoto and H. Kawamura, “Replica symmetry breaking in the RKKY skyrmion-crystal system,” Phys. Rev. B 104, 184432 (2021).
Mézard et al. (1987) M. Mézard, G. Parisi, and M. A. Virasoro, Spin glass theory and beyond: An Introduction to the Replica Method and Its Applications, Vol. 9 (World Scientific Publishing Company, 1987).
Mydosh (1993) J. A Mydosh, Spin glasses: an experimental introduction (CRC Press, 1993).
Kawamura and Taniguchi (2015) H. Kawamura and T. Taniguchi, “Spin glasses,” in Handbook of magnetic materials, Vol. 24 (Elsevier, 2015) pp. 1–137.
Ewald (1921) P. P. Ewald, “Die Berechnung optischer und elektrostatischer Gitterpotentiale,” (1921).
Hansen (1973) J. P. Hansen, “Statistical mechanics of dense ionized matter. I. Equilibrium properties of the classical one-component plasma,” Phys. Rev. A 8, 3096–3109 (1973).
Fuchizaki (1994) K. Fuchizaki, “Towards generalization of Ewald sum,” J. Phys. Soc. Jpn. 63, 4051–4059 (1994).
Ikeda and Kawamura (2008) A. Ikeda and H. Kawamura, “Ordering of the pyrochlore Ising model with the long-range RKKY interaction,” J. Phys. Soc. Jpn. 77, 073707–073707 (2008).
Creutz et al. (1979) M. Creutz, L. Jacobs, and C. Rebbi, “Monte Carlo study of Abelian lattice gauge theories,” Phys. Rev. D 20, 1915–1922 (1979).
Mermin and Wagner (1966) N. D. Mermin and H. Wagner, “Absence of ferromagnetism or antiferromagnetism in one- or two-dimensional isotropic Heisenberg models,” Phys. Rev. Lett. 17, 1133–1136 (1966).
Bruno (2001) P. Bruno, “Absence of spontaneous magnetic order at nonzero temperature in one- and two-dimensional Heisenberg and $\mathit{XY}$ systems with long-range interactions,” Phys. Rev. Lett. 87, 137203 (2001).
Binder, K. and Young, A. P. (1986) Binder, K. and Young, A. P., “Spin glasses: Experimental facts, theoretical concepts, and open questions,” Rev. Mod. Phys. 58, 801–976 (1986).
Parisi, G. and Urbani, P. and Zamponi, F. (2020) Parisi, G. and Urbani, P. and Zamponi, F., Theory of Simple Glasses: Exact Solutions in Infinite Dimensions (Cambridge University Press, 2020).
(38) Osamura, R. and Aoyama, K. and Mitsumoto, K. and Kawamura, H., in preparation.
Inosov et al. (2009) D. S. Inosov, D. V. Evtushinsky, A. Koitzsch, V. B. Zabolotnyy, S. V. Borisenko, A. A. Kordyuk, M. Frontzek, M. Loewenhaupt, W. Löser, I. Mazilu, H. Bitterlich, G. Behr, J.-U. Hoffmann, R. Follath, and B. Büchner, “Electronic structure and nesting-driven enhancement of the RKKY interaction at the magnetic ordering propagation vector in ${\mathrm{Gd}}_{2}{\mathrm{PdSi}}_{3}$ and ${\mathrm{Tb}}_{2}{\mathrm{PdSi}}_{3}$ ,” Phys. Rev. Lett. 102, 046401 (2009).
Fisher (1983) M. E. Fisher, “Scaling, universality and renormalization group theory,” in Critical Phenomena, edited by F. J. W. Hahne (Springer Berlin Heidelberg, Berlin, Heidelberg, 1983) pp. 1–139.
Kawamura (1988) H. Kawamura, “Renormalization-group analysis of chiral transitions,” Phys. Rev. B 38, 4916–4928 (1988).
Oike, H. and Kikkawa, A. and Kanazawa, N. and Taguchi, Y. and Kawasaki, M. and Tokura, Y. and Kagawa, F. (2016) Oike, H. and Kikkawa, A. and Kanazawa, N. and Taguchi, Y. and Kawasaki, M. and Tokura, Y. and Kagawa, F., “Interplay between topological and thermodynamic stability in a metastable magnetic skyrmion lattice,” Nat. Phys. 12, 62–66 (2016).
Karube, K. and White, J. S. and Reynolds, N. and Gavilano, J. L. and Oike, H. and Kikkawa, A. and Kagawa, F. and Tokunaga, Y. and Rønnow, H. M. and Tokura, Y. and Taguchi, Y. (2016) Karube, K. and White, J. S. and Reynolds, N. and Gavilano, J. L. and Oike, H. and Kikkawa, A. and Kagawa, F. and Tokunaga, Y. and Rønnow, H. M. and Tokura, Y. and Taguchi, Y., “Robust metastable skyrmions and their triangular–square lattice structural transition in a high-temperature chiral magnet,” Nat. Mater. 15, 1237–1242 (2016).
Marmorini and Momoi (2014) Giacomo Marmorini and Tsutomu Momoi, “Magnon condensation with finite degeneracy on the triangular lattice,” Phys. Rev. B 89, 134425 (2014).
Lohani, V. and Hickey, C. and Masell, J. and Rosch, A. (2019) Lohani, V. and Hickey, C. and Masell, J. and Rosch, A., “Quantum Skyrmions in Frustrated Ferromagnets,” Phys. Rev. X 9, 041063 (2019).
Mallik et al. (1998) R. Mallik, E. V. Sampathkumaran, P. L. Paulose, H. Sugawara, and H. Sato, “Magnetic anomalies in gd2pdsi3,” Pramana 51, 505–509 (1998).