Effective bilinear-biquadratic model for noncoplanar ordering in itinerant magnets

TL;DR

This paper develops an effective bilinear-biquadratic spin model to explain noncoplanar magnetic orderings in itinerant magnets, capturing complex multiple-Q states and their dependence on spin-charge coupling and external magnetic fields.

Contribution

It introduces a minimal effective spin model beyond RKKY theory that explains the emergence of noncoplanar multiple-Q states in itinerant magnets.

Findings

The model reproduces vortex crystal and skyrmion crystal states.

Multiple-Q states are stabilized even with weak spin-charge coupling.

External magnetic fields induce a variety of complex magnetic orders.

Abstract

Noncollinear and noncoplanar magnetic textures including skyrmions and vortices act as emergent electromagnetic fields and give rise to novel electronic and transport properties. We here report a unified understanding of noncoplanar magnetic orderings emergent from the spin-charge coupling in itinerant magnets. The mechanism has its roots in effective multiple spin interactions beyond the conventional Ruderman-Kittel-Kasuya-Yosida (RKKY) mechanism, which are ubiquitously generated in itinerant electron systems with local magnetic moments. Carefully examining the higher-order perturbations in terms of the spin-charge coupling, we construct a minimal effective spin model composed of the bilinear and biquadratic interactions with particular wave numbers dictated by the Fermi surface. We find that our effective model captures the underlying physics of the instability toward noncoplanar multiple-$Q$ states discovered recently in itinerant magnets: a single-$Q$ helical state expected from the RKKY theory is replaced by a double-$Q$ vortex crystal with chirality density waves even for an infinitely small spin-charge coupling on generic lattices [R. Ozawa $\mathit{et \ al.}$, J. Phys. Soc. Jpn. $\mathbf{85}$, 103703 (2016)], and a triple-$Q$ skyrmion crystal with a high topological number of two appears with increasing the spin-charge coupling on a triangular lattice [R. Ozawa, S. Hayami, and Y. Motome, to appear in Phys. Rev. Lett. (arXiv: 1703.03227)]. We find that, by introducing an external magnetic field, our effective model exhibits a plethora of multiple-$Q$ states. Our findings will serve as a guide for exploring further exotic magnetic orderings in itinerant magnets.

Figures

Equations

$$\mathcal{H}=-\sum _{i,j,\sigma }{t}_{ij}{c}_{i\sigma }^{†}{c}_{j\sigma }+J\sum _{i,\sigma ,{\sigma }^{\prime }}{c}_{i\sigma }^{†}{\mathbf{\sigma }}_{\sigma {\sigma }^{\prime }}{c}_{i{\sigma }^{\prime }}\cdot {\mathbf{S}}_i,$$

$$\mathcal{H}=\sum _{\mathbf{k},\sigma }{\epsilon }_{\mathbf{k}}{c}_{\mathbf{k}\sigma }^{†}{c}_{\mathbf{k}\sigma }+\frac{J}{\sqrt{N}}\sum _{\mathbf{k},\mathbf{q},\sigma ,{\sigma }^{\prime }}{c}_{\mathbf{k}\sigma }^{†}{\mathbf{\sigma }}_{\sigma {\sigma }^{\prime }}{c}_{\mathbf{k}+\mathbf{q}{\sigma }^{\prime }}\cdot {\mathbf{S}}_{\mathbf{q}},$$

$${\epsilon }_{\mathbf{k}}=-2\sum _{l=1,2}(t_1\cos \mathbf{k}\cdot {\mathbf{e}}_l+t_3\cos 2\mathbf{k}\cdot {\mathbf{e}}_l),$$

$${\epsilon }_{\mathbf{k}}=-2\sum _{l=1,2,3}(t_1\cos \mathbf{k}\cdot {\mathbf{e}}_l+t_3\cos 2\mathbf{k}\cdot {\mathbf{e}}_l),$$

$$F-F^{(0)}$$

$$=-\frac{T}{2!}{\int }_0^{\beta }d{\tau }_1{\int }_0^{\beta }d{\tau }_2⟨\mathcal{T}{\mathcal{H}}^{\prime }({\tau }_1){\mathcal{H}}^{\prime }({\tau }_2){⟩}_{\mathrm{con}}$$

$$-\frac{T}{4!}{\int }_0^{\beta }d{\tau }_1\cdots {\int }_0^{\beta }d{\tau }_4⟨\mathcal{T}{\mathcal{H}}^{\prime }({\tau }_1)\cdots {\mathcal{H}}^{\prime }({\tau }_4){⟩}_{\mathrm{con}}$$

$$-\cdots$$

$$=F^{(2)}+F^{(4)}+\cdots ,$$

$$F^{(2n)}$$

$$F^{(2)}=T\frac{J^2}{N}\sum _{\mathbf{k},\mathbf{q},{\omega }_p}{G}_{\mathbf{k}+\mathbf{q}}{G}_{\mathbf{k}}{\mathbf{S}}_{\mathbf{q}}\cdot {\mathbf{S}}_{-\mathbf{q}}.$$

$$F^{(2)}$$

$${\chi }_{\mathbf{q}}^0$$

$$=\frac{1}{N}\sum _{\mathbf{k}}\frac{f({\epsilon }_{\mathbf{k}})-f({\epsilon }_{\mathbf{k}+\mathbf{q}})}{{\epsilon }_{\mathbf{k}+\mathbf{q}}-{\epsilon }_{\mathbf{k}}}.$$

$${\mathbf{S}}_i=(\cos \mathbf{Q}\cdot {\mathbf{r}}_i,\sin \mathbf{Q}\cdot {\mathbf{r}}_i,0).$$

$$F^{(2)}$$

$$F^{(4)}=\frac{T}{2}\frac{J^4}{N^2}\sum _{\mathbf{k},{\omega }_p}\sum _{{\mathbf{q}}_1,{\mathbf{q}}_2,{\mathbf{q}}_3,{\mathbf{q}}_4,l}$$

$$\times \left[({\mathbf{S}}_{{\mathbf{q}}_1}\cdot {\mathbf{S}}_{{\mathbf{q}}_2})({\mathbf{S}}_{{\mathbf{q}}_3}\cdot {\mathbf{S}}_{{\mathbf{q}}_4})+({\mathbf{S}}_{{\mathbf{q}}_1}\cdot {\mathbf{S}}_{{\mathbf{q}}_4})({\mathbf{S}}_{{\mathbf{q}}_2}\cdot {\mathbf{S}}_{{\mathbf{q}}_3})-({\mathbf{S}}_{{\mathbf{q}}_1}\cdot {\mathbf{S}}_{{\mathbf{q}}_3})({\mathbf{S}}_{{\mathbf{q}}_2}\cdot {\mathbf{S}}_{{\mathbf{q}}_4})\right].$$

$$F^{(4)}$$

$$F_1^{(4)}$$

$$F_2^{(4)}$$

$$F_3^{(4)}$$

$$F_4^{(4)}$$

$$F_5^{(4)}$$

$$A_1$$

$$A_2$$

$$B_1$$

$$B_2$$

$$B_3$$

$$F_{(\mathbf{Q},-\mathbf{Q})}^{(2n)}=$$

$$\times \sum _{\mathbf{k},{\omega }_p,\nu }(G_{\mathbf{k}}{)}^n(G_{\mathbf{k}+{\mathbf{Q}}_{\nu }}{)}^n({\mathbf{S}}_{{\mathbf{Q}}_{\nu }}\cdot {\mathbf{S}}_{-{\mathbf{Q}}_{\nu }}{)}^n.$$

$$F_{(\mathbf{Q},-\mathbf{Q})}$$

$$=2T\frac{J^2}{N}\sum _{\mathbf{k},{\omega }_p,\nu }\left[\frac{G_{\mathbf{k}}{G}_{\mathbf{k}+{\mathbf{Q}}_{\nu }}}{1-2(J^2/N)G_{\mathbf{k}}{G}_{\mathbf{k}+{\mathbf{Q}}_{\nu }}({\mathbf{S}}_{{\mathbf{Q}}_{\nu }}\cdot {\mathbf{S}}_{-{\mathbf{Q}}_{\nu }})}\right]({\mathbf{S}}_{{\mathbf{Q}}_{\nu }}\cdot {\mathbf{S}}_{-{\mathbf{Q}}_{\nu }}).$$

$$\mathcal{H}=2\sum _{\nu }\left[-\tilde{J}{\mathbf{S}}_{{\mathbf{Q}}_{\nu }}\cdot {\mathbf{S}}_{-{\mathbf{Q}}_{\nu }}+\tilde{K}({\mathbf{S}}_{{\mathbf{Q}}_{\nu }}\cdot {\mathbf{S}}_{-{\mathbf{Q}}_{\nu }}{)}^2\right],$$

$$S_s(\mathbf{q})=S_s^{xx}(\mathbf{q})+S_s^{yy}(\mathbf{q})+S_s^{zz}(\mathbf{q}),$$

$$S_s^{\alpha \alpha }(\mathbf{q})$$

$$S_s^{\perp }(\mathbf{q})=S_s^{xx}(\mathbf{q})+S_s^{yy}(\mathbf{q}),$$

$$m_{\mathbf{q}}=\sqrt{S_s(\mathbf{q})/N}.$$

$$S_{{\chi }^{\mathrm{sc}}}(\mathbf{q})=\frac{1}{N}\sum _{i,j}⟨{\chi }_i^{\mathrm{sc}}{\chi }_j^{\mathrm{sc}}⟩e^{i\mathbf{q}\cdot ({\mathbf{r}}_i-{\mathbf{r}}_j)},$$

Full text

††preprint: APS/123-QED

I Introduction

In condensed matter physics, it is a central issue to explore unusual electronic orderings because they bring a major advance in the fundamental physics and open a new path toward applications to next-generation electronics and spintronics devices. Among them, magnetic orderings have been intensively investigated from both theory and experiment. Various fascinating behaviors have been observed for unusual magnetic structures, such as magnetoelectric effects in noncollinear magnets Katsura et al. (2005); Mostovoy (2006) and topological Hall effects in noncoplanar magnets Loss and Goldbart (1992); Ye et al. (1999); Ohgushi et al. (2000). From a theoretical point of view, the important issue is to clarify what kind of interaction is responsible for realizing such unusual magnetic orderings. A well-known example is the exchange interaction between localized magnetic moments in Mott insulators. Recent studies have revealed that such exchange interactions may result in unusual magnetic states, e.g., a skyrmion crystal with a noncoplanar topological spin texture Seki et al. (2012); Adams et al. (2012); Nagaosa and Tokura (2013), in the presence of the Dzyaloshinskii-Moriya (DM) interaction Dzyaloshinsky (1958); Moriya (1960) coming from the spin-orbit coupling Rößler et al. (2006); Yi et al. (2009); Binz et al. (2006) or frustration between competing interactions Okubo et al. (2012); Leonov and Mostovoy (2015); Lin and Hayami (2016); Hayami et al. (2016a); Batista et al. (2016).

Another representative magnetic interaction has been discussed for peculiar magnetism in rare-earth and other itinerant magnets where localized magnetic moments interact with itinerant electrons. The exchange coupling to itinerant electrons gives rise to an effective magnetic interaction between localized moments, which is called the Ruderman-Kittel-Kasuya-Yosida (RKKY) interaction Ruderman and Kittel (1954); Kasuya (1956); Yosida (1957). In contrast to the short-ranged interactions in Mott insulators, the RKKY interaction is long-ranged (power-law decay) with changing the sign depending on the distance. It often leads to an instability toward helical ordering with a long period structure featured by a single-$Q$ modulation.

On the other hand, recent theoretical studies have shown that a single-$Q$ helical ordering in itinerant magnets is not stable and taken over by multiple-$Q$ modulated structures Martin and Batista (2008); Akagi and Motome (2010); Hayami and Motome (2014); Ozawa et al. (2016a). Such an instability originates from higher-order spin interactions beyond the RKKY interactions, appearing after tracing out the degrees of freedom of itinerant electrons Akagi et al. (2012); Hayami and Motome (2014); Ozawa et al. (2016a); Hayami et al. (2016b). For example, a triple-$Q$ noncoplanar magnetic order showing the topological Hall effect is realized at particular electron fillings on a triangular lattice, owing to effective four-spin interactions Martin and Batista (2008); Akagi and Motome (2010). Similar mechanism stabilizes double-$Q$ noncoplanar magnetic orders accompanying chirality density waves at generic filling on generic lattices Ozawa et al. (2016a).

Above two examples illustrate the importance of effective multiple spin interactions emergent in itinerant magnets in realizing multiple-$Q$ noncoplanar states rather than single-$Q$ helical states. The interplay between charge and spin degrees of freedom, however, generates a variety of multiple spin interactions, whose analysis becomes more complicated in the higher-order terms. Thus, it is crucial to elucidate essential contributions in itinerant magnets and construct a concise effective model for further exploration of exotic magnetic orderings. Such an effective model will be helpful to avoid laborious calculations for the itinerant electron systems, which in general need a tremendous computational cost.

In the present study, we clarify the minimal ingredient to capture the essential physics of exotic magnetism brought by the spin-charge coupling, focusing on the weak coupling regime. We show that the origin of noncoplanar magnetic orderings in itinerant magnets is an effective biquadratic interaction specified by particular wave numbers dictated by the Fermi surface. We derive an effective spin model with bilinear and biquadratic couplings by examining the dominant contributions in the perturbative expansion in terms of the spin-charge coupling. By constructing the phase diagram of the effective spin model on square and triangular lattices by Monte Carlo simulations, we find that our model provides a unified understanding of unconventional multiple-$Q$ magnetic orders previously found in itinerant magnets Ozawa et al. (2016a, b). We confirm the stability of the multiple-$Q$ phases obtained in the spin model by comparing with those in the original itinerant electron model by variational calculations. We also elucidate the magnetic phase diagram by applying an external magnetic field to our effective model. We find a variety of field-induced multiple-$Q$ phases including different types of double(triple)-$Q$ states on the square (triangular) lattice.

The rest of the paper is organized as follows. In Sec. II, we present the starting itinerant electron model, the Kondo lattice model. We discuss how effective multiple spin interactions are generated from the perturbative expansion for the Kondo lattice Hamiltonian with respect to the exchange coupling between itinerant electrons and localized spins. By carefully examining contributions from the expansion, we extract a minimal ingredient relevant to exotic magnetic orderings. In Sec. III, we construct the effective spin model with bilinear-biquadratic interactions defined in momentum space. In Sec. IV, we discuss the multiple-$Q$ instability in this effective model. In Sec. V, we compare the results between the effective spin model and the original Kondo lattice model by using variational calculations. In Sec. VI, we show the phase diagram under an external magnetic field. A plethora of multiple-$Q$ phases is obtained from our Monte Carlo simulations. Section VII is devoted to a summary and a discussion of candidate materials.

II Effective Multiple Spin Interactions in Itinerant magnets

In this section, we discuss effective exchange interactions between localized spins based on the perturbative expansion with respect to the exchange coupling in the Kondo lattice model. In Sec. II.1, we introduce the Hamiltonian of the model. In Sec. II.2, we present a general expression for the perturbative expansion in the Kondo lattice model. Then, in Sec. II.3, we discuss the effect of the second-order RKKY interaction, which is not enough to determine magnetic orderings even when the exchange coupling is infinitely small. We extract the minimal ingredient to induce noncoplanar magnetic orderings by examining the fourth-order spin interactions in Sec. II.4, and generalize it to higher orders in Sec. II.5.

II.1 Model

We begin with a Kondo lattice model consisting of itinerant electrons and localized spins on the square and triangular lattices. The Hamiltonian is given by

[TABLE]

where $c^{\dagger}_{i\sigma}$ ($c_{i\sigma}$) is a creation (annihilation) operator of an itinerant electron at site $i$ and spin $\sigma$. The first term represents the kinetic motion of itinerant electrons. We consider hopping elements between nearest-neighbor sites, $t_{ij}=t_{1}$, and third-neighbor sites, $t_{ij}=t_{3}$, in the following analyses. It is noteworthy that qualitative features derived from the model in Eq. (1) are expected to hold for other choices of the hopping elements, e.g., second-neighbor hopping instead of $t_{3}$, whenever the bare magnetic susceptibility shows multiple maxima at symmetry-related wave numbers, as discussed in Sec. II.3. Hereafter, we take $t_{1}=1$ as an energy unit of the model in Eq. (1). The second term represents the exchange coupling between itinerant electron spins and localized spins. $\bm{\sigma}=(\sigma^{x},\sigma^{y},\sigma^{z})$ is the vector of Pauli matrices, $\mathbf{S}_{i}$ is a localized spin at site $i$ which is regarded as a classical spin with length $|\mathbf{S}_{i}|=1$, and $J$ is the exchange coupling constant; the sign of $J$ is irrelevant for the classical treatment of $\mathbf{S}_{i}$.

For the following arguments, it is useful to express the Hamiltonian in Eq. (1) in momentum space as

[TABLE]

where $c_{\mathbf{k}\sigma}^{\dagger}$ and $c_{\mathbf{k}\sigma}$ are the Fourier transform of $c_{i\sigma}^{\dagger}$ and $c_{i\sigma}$, respectively. $\varepsilon_{\mathbf{k}}$ is the energy dispersion of free electrons depending on the lattice structures: for the square lattice,

[TABLE]

where $\mathbf{e}_{1}=\hat{\mathbf{x}}=(1,0)$ and $\mathbf{e}_{2}=\hat{\mathbf{y}}=(0,1)$, and for the triangular lattice,

[TABLE]

where $\mathbf{e}_{1}=\hat{\mathbf{x}}$, $\mathbf{e}_{2}=-\hat{\mathbf{x}}/2+\sqrt{3}\hat{\mathbf{y}}/2$, and $\mathbf{e}_{3}=-\hat{\mathbf{x}}/2-\sqrt{3}\hat{\mathbf{y}}/2$. We set the lattice constant $a=1$ as the length unit. In the second term in Eq. (2), $\mathbf{S}_{\mathbf{q}}$ is the Fourier transform of $\mathbf{S}_{i}$ and $N$ is the number of sites. The second term can be regarded as the scattering of itinerant electrons by localized spins with the momentum transfer $\mathbf{q}$.

II.2 Perturbation expansion

We consider effective spin interactions mediated by itinerant electrons in the model in Eq. (2). Suppose the exchange coupling $J$ is small enough compared to the bandwidth of itinerant electrons, one can expand the free energy of the system with respect to $J$:

[TABLE]

where $\mathcal{H}^{\prime}$ represents the second term of Eq. (2), $\mathcal{T}$ is the time-ordered product, $\tau$ is the imaginary time, $T$ is the temperature, and $\beta$ is the inverse temperature ($=1/T$; we set the Boltzmann constant as unity). $\langle\cdots\rangle_{\rm con}$ stands for the averaged value over the connected Feynman diagrams. $F^{(0)}$ represents the free energy from the first term of Eq. (2). Note that there is no odd term in the expansion due to the time reversal symmetry of the model. The $2n$th-order contribution in the free energy can be expressed in the general form

[TABLE]

where $G_{\mathbf{k}}(i\omega_{p})=\left[i\omega_{p}-(\varepsilon_{\mathbf{k}}-\mu)\right]^{-1}$ is noninteracting Green’s function, $\omega_{p}$ is the Matsubara frequency, $\mu$ is the chemical potential, $\delta$ is the Kronecker delta, and $\mathbf{G}$ is the reciprocal lattice vector ($l$ is an integer). The sum of $\{P\}$ is taken for all the combinations of $\nu$ and $\nu^{\prime}$ (the number of the combinations is ${}_{2n}\mathrm{C}_{2}\cdot{}_{2n-2}\mathrm{C}_{2}\cdots{}_{2}\mathrm{C}_{2}/(n!)$), and $\lambda_{P}$ is $+1$ ($-1$) for even (odd) permutation. The product is taken for $1<\nu^{\prime}<\nu<2n$; see also the explicit expressions for the case with $n=1$ and $n=2$ in Eqs. (9) and (II.4), respectively. In Eq. (8), we omit the spin dependence of Green’s function because the unperturbed Hamiltonian is independent of the spin index. The similar expressions were discussed in Ref. Komarov and Dzebisashvili, 2016.

Figure 1 represents the Feynman diagrams for the first three terms in the expansion of Eq. (5), namely, $n=1$, $2$, and $3$ in Eq. (8). The diagrams consist of the scattering vertices by localized spins and the bare propagators of itinerant electrons, $G_{\mathbf{k}}$. By taking summations in terms of $\mathbf{q}_{\nu}$ ($\nu=1$, $2$, $\cdots$, $2n$) in Eq. (8), we can obtain the multiple spin interactions at any order of $J$. In the following sections, we discuss the specific form of such interactions by focusing on the second order (Sec. II.3) and fourth order (Sec. II.4) of $J$, and give the general expression for the higher-order contributions (Sec. II.5). Hereafter, we do not explicitly indicate the Matsubara frequency dependence in Green’s functions for simplicity.

II.3 Second-order RKKY interaction

First, let us consider the lowest-order contribution in Eq. (8), i.e., the second order of $J$ ($n=1$). This is obtained by substituting $n=1$ into Eq. (8), which is explicitly written as

[TABLE]

By taking the summation of $\omega_{p}$, Eq. (9) turns into

[TABLE]

where $\chi_{\mathbf{q}}^{0}$ is the bare susceptibility of itinerant electrons,

[TABLE]

Eq. (10) gives a pairwise interaction between localized spins, which is called the RKKY interaction Ruderman and Kittel (1954); Kasuya (1956); Yosida (1957). The sign and amplitude of the interaction depend on the band structure and electron filling through Eq. (12).

The magnetic state that optimizes the RKKY interaction in Eq. (10) is a single-$Q$ helical (spiral) state, whose spin structure is represented by

[TABLE]

Here, $\mathbf{Q}$ is the ordering vector defining the pitch and direction of the spiral, which is dictated by the peak of $\chi_{\mathbf{q}}^{0}$ in Eq. (12). Note that the spiral axis in the helical state in Eq. (13) is arbitrary because of the spin rotational symmetry of the RKKY interaction in Eq. (10). The reason why the helical state is preferred is understood from the normalization condition $|\mathbf{S}_{i}|=1$, which is equivalent to $\sum_{\mathbf{q}}|\mathbf{S}_{\mathbf{q}}|^{2}=N$: the helical state, in which $|\mathbf{S}_{\mathbf{Q}}|^{2}=|\mathbf{S}_{-\mathbf{Q}}|^{2}=N/2$ and $\mathbf{S}_{\mathbf{q}}=0$ for $\mathbf{q}\neq\pm\mathbf{Q}$, gives the lowest energy of Eq. (10). In other words, any superposition with another wave number or any higher harmonics leads to an energy cost compared to the helical state. Thus, the second-order free energy for the helical state is given by

[TABLE]

The corresponding Feynman diagram is shown in Fig. 2.

An important point in the RKKY analysis is that there are still four(six)fold degeneracy in the helical state on the square (triangular) lattice. Among them, the twofold degeneracy comes from the chiral symmetry in the centrosymmetric Bravais lattice structures, while the remaining two(three)fold degeneracy from the $C_{4}$ ($C_{6}$) rotational symmetry of the square (triangular) lattice structures. The latter degeneracy plays a role in realizing multiple-$Q$ orderings as discussed in Sec. II.4.

Let us describe how the degeneracy arises from the rotational symmetry by showing the momentum dependence of the bare susceptibility. $\chi_{\mathbf{q}}^{0}$ possesses multiple peaks reflecting the rotational symmetry of the lattice structure. This is demonstrated in Fig. 3: Figs. 3(a) and 3(c) show $\chi_{\mathbf{q}}^{0}$ on the square lattice with $t_{3}=-0.5$ and $\mu=0.98$ and the triangular lattice with $t_{3}=-0.85$ and $\mu=-3.5$, respectively. The corresponding Fermi surfaces are also shown in Figs. 3(b) and 3(d). The bare susceptibility shows multiple peaks at the wave numbers for which the Fermi surface is nested, and the wave numbers respect the rotational symmetry of the system. The former square lattice case has four peak structures at $\mathbf{Q}_{1}=\pm(2\pi/6,2\pi/6)$ and $\mathbf{Q}_{2}=R(\pi/2)\mathbf{Q}_{1}$, and the latter triangular lattice case shows six peak structures at $\mathbf{Q}_{1}=(2\pi/6,0)$, $\mathbf{Q}_{2}=R(2\pi/3)\mathbf{Q}_{1}$ and $\mathbf{Q}_{3}=R(4\pi/3)\mathbf{Q}_{1}$; here, $R(\theta)$ represents the rotational operator by $\theta$. As described above, at the level of the RKKY interaction in Eq. (10), the single-$Q$ helical order is realized by choosing the ordering vector $\mathbf{Q}$ out of these multiple peaks.

In this way, the second-order free energy, i.e., the RKKY interaction, favors the helical ordering with the single-$Q$ modulation represented by Eq. (13). However, the ground state in the Kondo lattice model is not given by the helical state but by noncoplanar multiple-$Q$ states even for the $J\to 0$ limit. The striking result was originally found at particular electron fillings where the Fermi surface has perfect nesting Martin and Batista (2008) or multiple connections Akagi and Motome (2010); Akagi et al. (2012); Hayami and Motome (2014), while recently generalized to generic fillings where the Fermi surface has no special property except for the rotational symmetry Ozawa et al. (2016a). The fundamental mechanism is that the system tends to lift the degeneracy due to the rotational symmetry of the lattice structure through the higher-order contributions of the free energy, as described in the following sections.

II.4 Fourth-order interaction

Next, we consider the fourth-order contribution of the free energy, Eq. (8) with $n=2$. It is explicitly written as

[TABLE]

This gives the four-spin interactions, which play an important role in lifting the degeneracy between the helical orderings and leads to the instability toward multiple-$Q$ orderings. For discussing such degeneracy lifting, it is enough to take into account the wave numbers for the multiple maxima in the bare susceptibility: $\mathbf{q}=\pm\mathbf{Q}_{1}$ and $\pm\mathbf{Q}_{2}$ ($\pm\mathbf{Q}_{1}$, $\pm\mathbf{Q}_{2}$, and $\pm\mathbf{Q}_{3}$) for the square (triangular) lattice. In the following, we consider the scattering processes satisfying $\mathbf{q}_{1}+\mathbf{q}_{2}+\mathbf{q}_{3}+\mathbf{q}_{4}=0$, i.e., $l=0$ in Eq. (II.4); the special cases with $\mathbf{q}_{1}+\mathbf{q}_{2}+\mathbf{q}_{3}+\mathbf{q}_{4}=\mathbf{G}$ were discussed for $2\mathbf{Q}_{\nu}=\mathbf{G}$ Martin and Batista (2008); Akagi and Motome (2010); Akagi et al. (2012); Hayami and Motome (2014, 2015) and $4\mathbf{Q}_{\nu}=\mathbf{G}$ Hayami et al. (2016b) ($\nu=1,2,3$). Then, the fourth-order free energy is given by the sum of five types of multiple spin interactions:

[TABLE]

where the sums in Eqs. (19)-(21) are taken for $\nu>\nu^{\prime}$. The coefficients are given by

[TABLE]

Each contribution is expressed by the diagram in Fig. 4. $A_{1}$ and $A_{2}$ represent the scattering processes by a single wave number, while $B_{1}$, $B_{2}$, and $B_{3}$ by multiple wave numbers.

Let us discuss which term plays a dominant role among the five types of multiple spin interactions in Eqs. (17)-(21). Figures 5(a) and 5(b) compares the coefficients $A_{1}$, $A_{2}$, $B_{1}$, $B_{2}$, and $B_{3}$ for two sets of parameters used in Fig. 3. In both cases, the coefficient $A_{2}$ becomes dominant in the low-temperature limit. This indicates that $F_{2}^{(4)}$ is the most important contribution among the fourth-order multiple spin interactions, as confirmed in Figs. 5(c) and 5(d).

It is worth noting that $F_{2}^{(4)}$ in Eq. (18) is in the form of the biquadratic interaction with a positive coefficient $A_{2}$. The positive biquadratic interaction, in general, favors a noncollinear spin configuration. For example, a triple-$Q$ noncoplanar ordering is realized on a triangular lattice, when the positive biquadratic interaction is enhanced by Fermi surface connections Akagi et al. (2012); Hayami and Motome (2014).

II.5 Higher-order contributions

The higher-order contributions can be straightforwardly expressed by the Feynman diagrams similar to those in Fig. 4. By analyzing the contributions from each diagram, we find that the scattering process proportional to $(\mathbf{S}_{\mathbf{Q}_{\nu}}\cdot\mathbf{S}_{-\mathbf{Q}_{\nu}})^{n}$ gives the most important contribution in each order. Note that the dominant $F_{2}^{(4)}$ in Eq. (18) is the case with $n=2$. Thus, the general form of the dominant contribution in the free energy at the $2n$th order is given by

[TABLE]

The sum of the dominant contributions up to the infinite orders can be summarized in a compact form

[TABLE]

The corresponding Feynman diagrams are shown in Fig. 6.

We find that the dominant terms in Eq. (28) contribute in a different way depending on the order of the expansion: the $(4l+2)$th-order terms proportional to $TG^{2l+1}_{\mathbf{k}}G^{2l+1}_{\mathbf{k}+\mathbf{Q}}<0$ tend to favor a single-$Q$ coplanar order, while the $4l$th-order ones proportional to $TG^{2l}_{\mathbf{k}}G^{2l}_{\mathbf{k}+\mathbf{Q}}>0$ tend to favor multiple-$Q$ noncoplanar order. This implies that Eq. (28) is divided into two groups, which are phenomenologically represented by a bilinear interaction $(\mathbf{S}_{\mathbf{Q}_{\nu}}\cdot\mathbf{S}_{-\mathbf{Q}_{\nu}})$ and a biquadratic interaction $(\mathbf{S}_{\mathbf{Q}_{\nu}}\cdot\mathbf{S}_{-\mathbf{Q}_{\nu}})^{2}$. On the basis of this observation, we construct an effective model in Sec. III.

We note that the free energy in Eq. (29) obtained the partial summation of the dominant contributions converges in the limit of zero temperature, while each term is divergent, as inferred in Figs. 5(c) and 5(d). We discuss the effect of the higher-order contributions by evaluating Eq. (29) in Sec. V.

III Effective spin model

III.1 Bilinear-biquadratic model in momentum space

The perturbation expansion in Sec. II indicates that many different types of effective spin interactions can contribute to the magnetic ordering in itinerant magnets. The careful comparison between different terms, however, gives an insight into the dominant contribution, as discussed in Secs. II.4 and II.5. Based on the observations, we propose an effective spin model in the weak coupling regime by including the contributions from the bilinear and biquadratic interactions. The Hamiltonian is given by

[TABLE]

where the sum is taken for $\nu=1,2$ ($1,2,3$) for the square (triangular) lattice, and $\mathbf{Q}_{\nu}$ are the wave numbers for the multiple peaks of the bare susceptibility $\chi_{\mathbf{q}}^{0}$  $\tilde{J}$ and $\tilde{K}$ are the coupling constants for bilinear and biquadratic interactions, respectively, in momentum space. We take $\tilde{J}>0$ and $\tilde{K}>0$, following the most dominant contributions to each term, Eqs. (14) and (18). Hereafter, we set $\tilde{J}=1$ and change the parameter $K=N\tilde{K}$ in the following calculations.

The model in Eq. (30) is a bilinear-biquadratic model defined in momentum space. The bilinear-biquadratic model has been studied in real space: the interactions are defined between the localized spins at particular sites, like $\mathbf{S}_{i}\cdot\mathbf{S}_{j}$ and $(\mathbf{S}_{i}\cdot\mathbf{S}_{j})^{2}$ in Refs. Sivardière and Blume, 1972; Chen and Levy, 1973; Harada and Kawashima, 2002. In our model, however, the interactions are defined for the Fourier components of spins whose wave numbers are dictated by the Fermi surface 111The effective model can be applicable to the case of $2\mathbf{Q}_{\nu}=\mathbf{G}$ and $4\mathbf{Q}_{\nu}=\mathbf{G}$ ($\nu=1,2,3$).

The perturbative argument in Sec. II suggests that the bilinear $\tilde{J}$ term is dominant over the biquadratic $\tilde{K}$ term, as the former originates from the RKKY interaction proportional to $J^{2}$ [Eq. (14)], while the most dominant contribution to the latter is proportional to $J^{4}$ [Eq. (18)] in the weak coupling limit. Nevertheless, we will study the effective model in Eq. (30) up to $K\sim\tilde{J}$ in the following sections, since the dominant contributions from the higher-order perturbation are phenomenologically renormalized into the effective bilinear and biquadratic interactions, as discussed in Sec. II.5. We will demonstrate that the extension of the model to the nonperturbative region $K\sim\tilde{J}$ is indeed useful to discuss the magnetic instabilities appearing in the original Kondo lattice model in Eq. (1).

III.2 Monte Carlo simulation

In order to obtain the magnetic phase diagram of the effective spin model in Eq. (30), we perform classical Monte Carlo simulations at low temperatures. Our simulations are carried out with the standard Metropolis local updates. In the following, we present the results for the systems with $N=96^{2}$ sites under periodic boundary conditions in both the square and triangular lattice cases; we confirm that the system size dependence is small by performing the simulations for $N=48^{2}$ and $72^{2}$. In each simulation, we perform simulated annealing to find the low-energy configuration in the following way. We gradually reduce the temperature with a rate $T_{n+1}=\alpha T_{n}$ where $T_{n}$ is the temperature in the $n$th step. We set the initial temperature $T_{0}=0.1$-$1.0$ and take the coefficient of geometrical cooling $\alpha=0.9995$-$0.9999$. The final temperature, which is typically taken at $T=0.01$, is reached by spending totally $10^{5}$-$10^{6}$ Monte Carlo sweeps. At the final temperature, we perform $10^{5}$-$10^{6}$ Monte Carlo sweeps for measurements after $10^{5}$-$10^{6}$ steps for thermalization.

We calculate the structure factors for spin and scalar chirality to identify each magnetic phase. The spin structure factor $S_{s}(\mathbf{q})$ is given by

[TABLE]

where

[TABLE]

with $\alpha=x,y,z$. We also compute

[TABLE]

in the presence of the external magnetic field applied along the $z$ direction. We also introduce the following notation for the magnetic moments at $\mathbf{q}$ components:

[TABLE]

In addition, we calculate the chirality structure factor. For the square lattice, it is defined by

[TABLE]

where the local scalar chirality at site $i$ is introduced as $\chi^{\rm sc}_{i}=\mathbf{S}_{i}\cdot(\mathbf{S}_{i+\hat{x}}\times\mathbf{S}_{i+\hat{y}})+\mathbf{S}_{i}\cdot(\mathbf{S}_{i-\hat{x}}\times\mathbf{S}_{i-\hat{y}})-\mathbf{S}_{i}\cdot(\mathbf{S}_{i-\hat{x}}\times\mathbf{S}_{i+\hat{y}})-\mathbf{S}_{i}\cdot(\mathbf{S}_{i+\hat{x}}\times\mathbf{S}_{i-\hat{y}})$. Meanwhile, the chirality structure factor on the triangular lattice is defined separately for the upward and downward triangles as

[TABLE]

where $\mathbf{R}$ and $\mathbf{R}^{\prime}$ represent the position vectors at the centers of triangles, and $\mu=(u,d)$ represent upward and downward triangles, respectively. Note that connecting the center of gravity in the upward and downward triangles leads to the honeycomb networks. In Eq. (36), $\chi^{\rm sc}_{\mathbf{R}}=\mathbf{S}_{j}\cdot(\mathbf{S}_{k}\times\mathbf{S}_{l})$, where $j,k,l$ are the sites on the triangle at $\mathbf{R}$ in the counterclockwise order. We also introduce the scalar chirality at $\mathbf{q}$ components:

[TABLE]

IV Multiple-$Q$ magnetic instability

In this section, we investigate magnetic instabilities appearing in the bilinear-biquadratic model in Eq. (30), which is derived as an effective model for the Kondo lattice model with small $J/t$. In Sec. IV.1, we elucidate the phase diagrams on the square and triangular lattices, in which the helical state expected from the RKKY interaction is no longer stable and replaced by multiple-$Q$ states. In Sec. IV.2, we discuss why the multiple-$Q$ states have lower energy than the single-$Q$ helical state.

IV.1 Phase diagram

IV.1.1 Square lattice

First, we show the result on the square lattice. In order to investigate the magnetic instability, we perform Monte Carlo simulations at sufficiently low temperature, $T=0.01$, with using the simulated annealing described in the previous section. We set the ordering vectors as $\mathbf{Q}_{1}=(2\pi/6,2\pi/6)$ and $\mathbf{Q}_{2}=R(\pi/2)\mathbf{Q}_{1}$ in the model in Eq. (30) defined on the square lattice, although the following results are qualitatively unchanged for other sets of symmetry-related $\mathbf{Q}_{\nu}$.

Figure 7(a) shows $K$ dependences of the $\mathbf{Q}_{\nu}$ component of the magnetization, $m_{\mathbf{Q}_{\nu}}$. At $K=0$, the single-$Q$ helical state is realized because Eq. (30) is reduced to the RKKY interaction in Eq. (10). Indeed, our Monte Carlo result indicates that $m_{\mathbf{Q}_{1}}$ becomes nonzero, while $m_{\mathbf{Q}_{2}}$ vanishes (the state with $m_{\mathbf{Q}_{1}}=0$ and $m_{\mathbf{Q}_{2}}\neq 0$ is also obtained depending on the initial conditions and the scheduling in the simulated annealing). Once we turn on $K$, the $\mathbf{Q}_{2}$ component becomes nonzero and develops as increasing $K$. This indicates that the single-$Q$ state has an instability toward the double-$Q$ state for an infinitesimal $K$. With increasing $K$, $m_{\mathbf{Q}_{1}}$ decreases and $m_{\mathbf{Q}_{2}}$ increases, and they gradually approaches the same value, while they take different values, at least, in the calculated range for $K<1.2$.

We present a Monte Carlo snapshot of the spin configuration of this double-$Q$ state in the left panel of Fig. 8. The spin configuration is clearly different from the single-$Q$ helical state and modulated by the second $Q$ component. The spin structure factor in the middle panel of Fig. 8 shows two pairs of Bragg peaks at $\pm\mathbf{Q}_{1}$ and $\pm\mathbf{Q}_{2}$ with different intensities.

On the other hand, in terms of the spin scalar chirality, the double-$Q$ state is accompanied by the chiral density wave only with the wave number $\mathbf{Q}_{2}$, as shown in Fig. 7(b). Indeed, the structure factor shown in the right panel of Fig. 8 has a peak only at $\pm\mathbf{Q}_{2}$.

We summarize the phase diagram in Fig. 7(c). As mentioned above, the single-$Q$ helical state is realized only at $K=0$, and the double-$Q$ state with chirality stripe is stable for all $K>0$: the critical value of $K$ is $K_{c}=0$. We call this double-$Q$ state for $K>0$ “double-$Q^{\prime}$ chiral stripe”, where $Q^{\prime}$ means that the amplitudes of the $\mathbf{Q}_{1}$ and $\mathbf{Q}_{2}$ components are different.

We find that this double-$Q^{\prime}$ chiral stripe is similar to the double-$Q$ state recently found in the Kondo lattice model in Ref. Ozawa et al., 2016a. The real-space spin structure was described in the form

[TABLE]

where $b$ represents the amplitude of the $\mathbf{Q}_{2}$ component; the superscript T denotes the transpose of the vector. As discussed in the previous study Ozawa et al. (2016a), the spin configuration in Eq. (41) has two notable aspects. One is that it is continuously connected to the single-$Q$ helical state as $b\to 0$, which indicates that the double-$Q$ state is topologically trivial. The other is that the $\mathbf{Q}_{2}$ component is introduced in the perpendicular direction to the $\mathbf{Q}_{1}$ helical plane. Due to the latter, the energy loss from the higher harmonics at the RKKY level is small enough compared to the energy gain in the higher-order terms Ozawa et al. (2016a). We will discuss this point for our effective model in Sec. IV.2.

IV.1.2 Triangular lattice

Next, we study the phase diagram on the triangular lattice. We consider the specific set of ordering vectors $\mathbf{Q}_{1}=(2\pi/6,0)$, $\mathbf{Q}_{2}=R(2\pi/3)\mathbf{Q}_{1}$, and $\mathbf{Q}_{3}=R(4\pi/3)\mathbf{Q}_{1}$ in the model in Eq. (30) defined on the triangular lattice, while qualitative features are expected to hold for arbitrary sets of symmetry-related $\mathbf{Q}_{\nu}$, as in the square lattice case. The obtained results of the magnetic moments and chirality are shown in Figs. 9(a) and 9(b), respectively. In the spin sector, the single-$Q$ state at $K=0$ turns into a double-$Q$ state for $K>0$. We find that this double-$Q$ state corresponds to a straightforward extension of the double-$Q^{\prime}$ chiral stripe obtained for the square lattice: the spin structure is expressed by the same form as Eq. (41) with choosing two $\mathbf{Q}_{\nu}$ out of three [e.g., $\mathbf{Q}_{1}$ and $\mathbf{Q}_{2}$ in the solution obtained in Figs. 9(a) and 9(b)].

While increasing $K$, however, we find another phase transition at $K\sim 0.37$. At the transition, the magnetic moment changes discontinuously, as shown in Fig. 9(a); all $m_{\mathbf{Q}_{\nu}}$ become nonzero for $K\gtrsim 0.37$, which indicates that the spin state is a triple-$Q$ order. There, the amplitudes of the $\mathbf{Q}_{1}$, $\mathbf{Q}_{2}$, and $\mathbf{Q}_{3}$ components are equivalent, as shown in Fig. 9(a). This is also clearly seen in the spin structure factor shown in Fig. 10. On the other hand, in the scalar chirality sector, the uniform ($\mathbf{q}=0$) component is induced, while the $\mathbf{Q}_{\nu}$ components all vanish, as shown in Figs. 9(b).

Thus, the state for $K\gtrsim 0.37$ is the triple-$Q$ order with nonzero net scalar chirality, which is common to the skyrmion states studied in noncentrosymmetric chiral magnets Rößler et al. (2006); Yi et al. (2009); Binz et al. (2006); Nagaosa and Tokura (2013) and centrosymmetric frustrated magnets Okubo et al. (2012); Leonov and Mostovoy (2015); Lin and Hayami (2016); Hayami et al. (2016a). However, it is qualitatively distinct from the ordinary skyrmion crystals. Actually, we find that the triple-$Q$ state is the same as the unconventional skyrmion crystal discovered in Ref. Ozawa et al., 2016b. The real-space spin structure is described by an equal superposition of three helices as

[TABLE]

where $N_{i}$ is the site-dependent renormalization factor to satisfy $|\mathbf{S}_{i}|=1$: $N_{i}=\sqrt{\sum_{\nu=1}^{3}\cos^{2}\mathbf{Q}_{\nu}\cdot\mathbf{r}_{i}}$. This triple-$Q$ state has mainly three different aspects from the conventional skyrmion crystals. The first point is that the topological number $n_{\rm sk}$ is two per magnetic unit cell, as shown in Fig. 10 Ozawa et al. (2016b). This is in sharp contrast to the previous skyrmions with $n_{\rm sk}=1$. The second point is the symmetry in the ordered phase. The obtained triple-$Q$ state preserves inversion and spin rotational symmetries [O(3) symmetry], while most of the previous skyrmions in centrosymmetric systems have only U(1) symmetry because they are stabilized in an applied magnetic field. This peculiar symmetry may lead to new Goldstone modes. The third is the absence of the nonzero $\mathbf{Q}$ components in the scalar chirality, as shown in Fig. 9(b). This is in contrast to the presence of the six peaks in the chirality structure factor in the skyrmion crystal phase, e.g., found in Ref. Hayami et al., 2016a. Following Ref. Ozawa et al., 2016b, we call this triple-$Q$ state the triple-$Q$ $n_{\rm sk}=2$ skyrmion crystal.

We summarize the phase diagram in Fig. 9(c). As in the square lattice case, the phase boundary between the single-$Q$ helical state and the double-$Q^{\prime}$ chiral stripe is at $K=K_{c1}=0$. In the triangular lattice case, however, there is another phase transition at $K=K_{c2}\sim 0.37$, above which the triple-$Q$ $n_{\rm sk}=2$ skyrmion crystal appears.

IV.2 Energy comparison

We here discuss why the multiple-$Q$ orderings in Eqs. (41) and (45) are realized in the presence of the biquadratic $K$ term, which mimics the fourth-order contribution in the Kondo lattice model. We discuss the energy comparison between the single-$Q$ helical state and the double-$Q^{\prime}$ chiral stripe in Sec. IV.2.1 and among the single-$Q$ helical state, the double-$Q^{\prime}$ chiral stripe, and triple-$Q$ $n_{\rm sk}=2$ skyrmion crystal in Sec. IV.2.2. We also comment on a coplanar double-$Q$ state on the square lattice, which is the counterpart of the triple-$Q$ $n_{\rm sk}=2$ skyrmion crystal on the triangular lattice, in Sec. IV.2.3.

IV.2.1 Double-$Q^{\prime}$ chiral stripe

First, we discuss the competition between the single-$Q$ helical and the double-$Q^{\prime}$ chiral stripe states. We here show that the single-$Q$ helical state has higher energy than the double-$Q^{\prime}$ chiral stripe for $0<K\ll\tilde{J}$ in our effective spin model in Eq. (30). The following argument is generic and applicable to both the square and triangular lattices. The argument is basically parallel to that for the original Kondo lattice model Ozawa et al. (2016a), which suggests that our effective spin model captures the instability toward the double-$Q^{\prime}$ chiral stripe in the Kondo lattice model.

In order to estimate the energy of the double-$Q^{\prime}$ chiral stripe, we expand the square root in Eq. (41) with respect to $b^{2}$:

[TABLE]

where $C_{2l}$ ($l=0,1,2,\cdots$) are given by

[TABLE]

and so on. Then, the spin pattern in Eq. (41) is represented by

[TABLE]

By Fourier transformation, we obtain

[TABLE]

By substituting Eqs. (51) and (52) into the effective spin model in Eq. (30), we obtain the energy per site for the double-$Q^{\prime}$ chiral stripe:

[TABLE]

In the second line, we neglect the contributions $Kb^{4}$ because we assume $\tilde{J}\gg K$.

The important observation in Eq. (56) is that there are no contributions proportional to $\tilde{J}b^{2}$. This comes from the fact that the modulation with the $\mathbf{Q}_{2}$ component is perpendicular to the helical plane with the $\mathbf{Q}_{1}$ component, as described in Sec. IV.1. Obviously, the second term in Eq. (56) leads to the energy loss in the double-$Q^{\prime}$ chiral stripe, which is consistent with the argument that the RKKY interaction favors the single-$Q$ helical state. On the other hand, the third term in Eq. (56) represents the energy gain in the double-$Q^{\prime}$ chiral stripe. The condition to stabilize the double-$Q^{\prime}$ chiral stripe, i.e., $E^{2Q^{\prime}}<E^{1Q}$ reads

[TABLE]

This means that an infinitesimal $K$ leads to the instability of the single-$Q$ helical state by introducing the second component with the amplitude $b$.

We note that, strictly speaking, we need to take into account the contributions from Eqs. (53) and (54) because they also give the energy in the order of $\tilde{J}b^{4}$. However, we can neglect such contributions when the bare susceptibility has considerably small amplitudes at $\mathbf{Q}_{1}\pm 2\mathbf{Q}_{2}$ compared to the distinct peaks at $\mathbf{Q}_{1}$ and $\mathbf{Q}_{2}$, as shown in Fig. 3. The detailed description including the contributions of higher harmonics was given in Ref. Ozawa et al., 2016a.

IV.2.2 Triple-$Q$ skyrmion crystal

Next, we turn to the stability of the triple-$Q$ state on the triangular lattice, which appears for relatively large $K$. We can rewrite the effective Hamiltonian in Eq. (30) into the form

[TABLE]

where $|S_{\mathbf{Q}_{\nu}}|^{2}=\mathbf{S}_{\mathbf{Q}_{\nu}}\cdot\mathbf{S}_{-\mathbf{Q}_{\nu}}$. From this form, for optimizing the energy with respect to $K=N\tilde{K}$, it is necessary to minimize the contribution from the second term in Eq. (IV.2.2) and to maximize that from the third term. This is accomplished by the conditions: (i) the amplitude of the triple-$Q$ component is taken to be the same because of the sum rule $\sum_{\nu}|S_{\mathbf{Q}_{\nu}}|^{2}=N$, and (ii) the directions of $\mathbf{Q}_{\nu}$ modulations are perpendicular to each other, namely, the magnetic moments for the $\mathbf{Q}_{\nu}$ components are perpendicular to one another. The latter condition indicates that the contribution of the higher harmonics is small, as discussed above for the double-$Q^{\prime}$ state in the square lattice case. The spin configuration in Eq. (45) satisfies these two conditions, and hence, it is chosen as the ground state in the large $K$ region.

Let us evaluate the energy of the triple-$Q$ $n_{\rm sk}=2$ skyrmion crystal. Although the comparison with the double-$Q^{\prime}$ chiral stripe is difficult because Eq. (56) is justified only for small $b$, we here compare the energy of the triple-$Q$ state with that of the single-$Q$ state. In the triple-$Q$ state in Eq. (45), the magnetic moments for the $\mathbf{Q}_{\nu}$ components are given by

[TABLE]

where $\eta$ represents a small correction from the higher harmonics, e.g., $\mathbf{Q}_{1}+2\mathbf{Q}_{2}$ and $3\mathbf{Q}_{1}$, leading to the energy loss at the level of the RKKY interaction. By substituting Eq. (59) into Eq. (30), we obtain the energy per site for the triple-$Q$ $n_{\rm sk}=2$ skyrmion crystal as

[TABLE]

Monte Carlo simulations give $\eta\simeq 0.05$ almost independent of $K$. From Eq. (60), we find that the energy of the triple-$Q$ $n_{\rm sk}=2$ skyrmion crystal is higher than the single-$Q$ state for small $K$. While increasing $K$, however, the energy of the triple-$Q$ $n_{\rm sk}=2$ skyrmion crystal becomes lower for $K>3/11\tilde{J}\sim 0.27\tilde{J}$ when we set $\eta=0.05$. This is qualitatively consistent with the result by numerical simulations in Fig. 9(a), giving the phase boundary at $K_{c2}\sim 0.37\tilde{J}$. The underestimate of $K_{c2}$ is reasonable as we neglect the double-$Q^{\prime}$ chiral stripe in the current analysis.

IV.2.3 Coplanar double-$Q$ state on the square lattice

Let us comment on the reason why the second transition at $K_{c2}$ is absent in the square lattice case. We can think of a counterpart of the triple-$Q$ state in Eq. (45) for the square lattice case, which may optimize the $K$ term along with the conditions discussed in Sec. IV.2.2. As there are two $\mathbf{Q}_{\nu}$ in the square lattice case, such an optimal state will be given by

[TABLE]

where $N_{i}=\sqrt{\sum_{\nu=1}^{2}\cos^{2}\mathbf{Q}_{\nu}\cdot\mathbf{r}_{i}}$. Although the coplanar double-$Q$ state is expected to be stabilized for $K\gg\tilde{J}$ on the square lattice as it satisfies the two conditions in Sec. IV.2.2, it is not obtained as the lowest-energy state in our Monte Carlo simulations.

Nonetheless, the coplanar double-$Q$ state sometimes appears as a metastable state when we perform simulated annealing with a sudden quench from high-temperature limit. We find that, although the energy of this coplanar double-$Q$ state is always higher than that in the double-$Q^{\prime}$ chiral stripe, the difference becomes smaller for larger $K$; they are almost degenerate for $K\gg\tilde{J}$ within the numerical accuracy. This is because the energy loss in the $K$ term for the coplanar double-$Q$ state is almost similar to that for the double-$Q^{\prime}$ chiral stripe when taking $b\sim 1$. The situation is in contrast to the triangular lattice case. The energy loss in the $K$ term for the triple-$Q$ $n_{\rm sk}=2$ skyrmion crystal is given by $(K/6)(1-4\eta)$ in Eq. (60), which is smaller than $(K/4)(1-4\eta)$ for the double-$Q^{\prime}$ chiral stripe and coplanar double-$Q$ states. Moreover, from the numerical simulations, we find that the effect of higher harmonics in the coplanar double-$Q$ state, mainly from $\pm\mathbf{Q}_{1}\pm 2\mathbf{Q}_{2}$ and $\pm 3\mathbf{Q}_{1}$, becomes larger compared to that in the double-$Q^{\prime}$ chiral stripe. Thus, the coplanar double-$Q$ state is not stabilized in the effective model in Eq. (30) on the square lattice, and hence, there is no second transition for this case.

V Comparison to Itinerant model

So far, we have investigated the phase diagram of the spin model in Eq. (30), which we propose as the effective model for the Kondo lattice model in Eq. (1) in the small $J$ region. In this section, we examine the validity of the effective model by comparing the obtained phase diagrams with those by the direct calculations for the Kondo lattice model.

V.1 Variational calculations

In order to confirm the validity of the effective spin model obtained from the perturbative argument, we examine the ground state of the original Kondo lattice model in Eq. (1). We here perform a variational calculation as follows: we compare the grand potential at zero temperature $\Omega=E-\mu n_{\rm e}$ for variational states with different magnetic orders in the localized spins, where $E=\langle\mathcal{H}\rangle/N$ is the internal energy per site, and $n_{\rm e}=\sum_{i\sigma}\langle c_{i\sigma}^{\dagger}c_{i\sigma}\rangle/N$ is the electron density. For the variational states, we assume the spin configurations given in Eqs. (13), (41), and (45) with the same wave numbers used in the effective spin model. In the double-$Q^{\prime}$ state in Eq. (41), we deal with $b$ as the variational parameter. For comparison, we also calculate the grand potential for the ferromagnetic state, which is expected for large $J$ by the double-exchange mechanism Zener (1951); Anderson and Hasegawa (1955); the spin configuration is given by $\mathbf{S}_{i}=(0,0,1)$. We consider these variational states in the system with $N=480^{2}$ under the periodic boundary conditions.

Figure 11(a) shows $J$ dependences of the grand potential at $t_{3}=-0.5$ and $\mu=0.98$ on the square lattice for the ferromagnetic state and the double-$Q^{\prime}$ chiral stripe measured from that for the single-$Q$ helical state. The similar calculations were done in Ref. Ozawa et al., 2016a. The results are consistent with those obtained in the effective spin model in Sec. IV. The grand potential for the double-$Q^{\prime}$ chiral stripe state is lower than that for the single-$Q$ helical state as well as the ferromagnetic state for all $J>0$ [see also the inset of Fig. 11(a) for small $J$].

The results for the triangular lattice are also consistent with those for the effective spin model. Figure 11(b) shows $J$ dependences of the grand potential at $t_{3}=-0.85$ and $\mu=-3.5$ measured from that for the single-$Q$ helical state. The double-$Q^{\prime}$ chiral stripe becomes the lowest-energy state for $0<J\lesssim 0.11$, while the triple-$Q$ $n_{\rm sk}=2$ skyrmion crystal state replaces it for $J\gtrsim 0.11$.

Thus, for both square and triangular lattice cases, the variational calculations indicate that the multiple-$Q$ states are stable for nonzero $J$ in the Kondo lattice model. The sequence of the phases is consistent with that in the effective spin model in Eq. (30) shown in Figs. 7(c) and 9(c). The results strongly support that the effective model captures the instabilities toward the multiple-$Q$ states inherent to the Kondo lattice model.

V.2 Higher-order corrections

Our effective spin model in Eq. (30) includes the two contributions inferred from the perturbative expansion of the Kondo lattice model in terms of $J$. One is the bilinear $\tilde{J}$ term representing the expansion terms favoring a single-$Q$ helical order, whose dominant contribution is from the second-order RKKY interaction in Eq. (9). The other is the biquadratic $\tilde{K}$ term representing the expansion terms favoring a multiple-$Q$ noncoplanar order, whose most dominant contribution is in the fourth-order terms, Eq. (18). On the other hand, as discussed in Sec. II.5, each term in the dominant contributions in the expansion, $(\mathbf{S}_{\mathbf{Q}_{\nu}}\cdot\mathbf{S}_{-\mathbf{Q}_{\nu}})^{n}$ [Eq. (II.5) and Fig. 6] diverges in the low-temperature limit, even though the sum up to infinite order in Eq. (29) converges. To reinforce the validity of our effective model formally including only the terms with $n=1$ and $2$, we here show that the infinite sum in Eq. (29) does not change the qualitative phase diagram on the square and triangular lattices obtained in Sec. IV.

Figure 12 shows $J$ dependences of $-F_{(\mathbf{Q},-\mathbf{Q})}/(2J^{2})$ in Eq. (29) for the single-$Q$ helical and the double-$Q^{\prime}$ chiral stripe on the square lattice [Fig. 12(a)] and for the single-$Q$ helical, the double-$Q^{\prime}$ chiral stripe, and the triple-$Q$ $n_{\rm sk}=2$ skyrmion crystal on the triangular lattice [Fig. 12(b)]. The parameters are taken to be the same as Fig. 11. For comparison, we also plot the $-(\Omega-\Omega^{(0)})/(2J^{2})$ obtained by the variational calculations, where $\Omega^{(0)}$ represents the grand potential for $J=0$. As shown in Fig. 12, the free energy including the higher-order dominant contributions also reproduces well the phase diagrams on both square and triangular lattices. The agreement supports the validity of our effective model including only the bilinear and biquadratic interactions in momentum space.

VI Effect of magnetic field

We have found that the effective spin model in Eq. (30) captures the magnetic instability in the Kondo lattice model in Eq. (1). The effective model is useful for exploration of further exotic magnetic orderings, as its simplicity allows us to bypass laborious calculations necessary for itinerant electron models. In this section, we extend our study to construct the magnetic phase diagram by incorporating the effect of an external magnetic field on the effective spin model. For this purpose, adding the Zeeman coupling term to the Hamiltonian in Eq. (30), we consider the Hamiltonian given by

[TABLE]

We discuss the magnetic phase diagrams under the magnetic field on the square lattice in Sec. VI.1 and the triangular lattice in Sec. VI.2.

VI.1 Square lattice

Figure 13 shows the $K$-$H$ phase diagram for the square lattice case obtained by our Monte Carlo simulations. We find four different phases other than the fully-polarized state in a large magnetic field. The phase boundaries are determined by analyzing the structure factors from Eqs. (31) to (36). All the phases show no net scalar chirality. In what follows, we describe the details of each phase one by one. We also discuss the magnetization processes.

double-$Q^{\prime}$ chiral stripe I (2$Q^{\prime}$-I)

This phase occupies a region at low $H$ and $K$. The spin configuration is characterized by the double-$Q$ noncoplanar modulation with different intensities, similar to the double-$Q^{\prime}$ chiral stripe at $H=0$ discussed in Fig. 8, but the magnetic field breaks the O(3) symmetry in the present state; see the typical spin configuration in Fig. 14(a). The $xy$ component of the spin structure factor has a single-$Q$ peak at $\mathbf{Q}_{1}$, while the $z$ component exhibits two peaks at $\mathbf{Q}_{1}$ and $\mathbf{Q}_{2}$ with different intensities, as shown in Fig. 14(a). The chirality structure factor shows a single-$Q$ peak at $\mathbf{Q}_{2}$, similar to the double-$Q^{\prime}$ chiral stripe at $H=0$.

double-$Q^{\prime}$ chiral stripe II (2$Q^{\prime}$-II)

This state occupies the low-$H$ and high-$K$ region of the phase diagram in Fig. 13 and appears next to the double-$Q^{\prime}$ chiral stripe I phase with increasing $K$. The spin configuration also resembles a double-$Q^{\prime}$ chiral stripe state at $H=0$ in Fig. 8, whereas not only $z$ but also $xy$ components of the spin structure factor are characterized by the double-$Q$ peak structure in the present case, as shown in Fig. 14(b). On the other hand, the chirality structure factor shows a single-$Q$ peak at $\mathbf{Q}_{2}$, similar to the chiral stripe at $H=0$. The transition from the double-$Q^{\prime}$ I to II is presumably caused by the tendency that larger $K$ favors a larger amplitude of the double-$Q$ components, as already seen in the absence of $H$, as shown in Fig. 7(a).

single-$Q$ conical state (1$Q$)

This state occupies a region at higher $H$ and lower $K$, appearing next to the double-$Q^{\prime}$ chiral stripe I phase with increasing $H$. This is a single-$Q$ conical state, whose spin configuration is given by $\mathbf{S}_{i}=(\sin\theta\cos\mathbf{Q}_{\nu}\cdot\mathbf{r}_{i},\sin\theta\sin\mathbf{Q}_{\nu}\cdot\mathbf{r}_{i},\cos\theta)$ with $\theta=\cos^{-1}(H/H_{\rm sat})$ where $H_{\rm sat}=2\tilde{J}$, independent of the value of $K$ [see Fig. 14(c)]. We note that the conical state is realized even at $K=0$, where the model is reduced to a frustrated Heisenberg model in a magnetic field. With increasing $H$, the conical angle $\theta$ becomes smaller, and this phase continuously changes into the fully-polarized state at $H=2$, as shown in Fig. 13.

double-$Q$ state (2$Q$)

This phase occupies the largest portion of the phase diagram in Fig. 13, which is found next to the single-$Q$ conical state upon increasing $K$. The spin configuration is similar to the double-$Q$ state in Eq. (64), although the present state shows a noncoplanar spin structure by spin canting along the magnetic field direction. It consists of a periodic array of vortices and antivortices, as shown in Fig. 14(d). The $xy$ component of the spin structure factor shows the double-$Q$ peaks with equal intensities, while the $z$ components are negligibly small, except for a large peak at $\mathbf{q}=0$ from the net moment induced by the magnetic field [omitted in Fig. 14(d) for clarity]. With increasing $H$, the phase continuously changes into the fully-polarized state at $H=2$, as the single-$Q$ conical state. Note that this double-$Q$ state is also similar to the one reported on the triangular lattice in classical Hayami et al. (2016c) and quantum Kamiya and Batista (2014); Marmorini and Momoi (2014) magnets.

Magnetization curves

Figure 15 shows the magnetization curves at several $K$. The magnetization is given by $M=m_{\mathbf{q}=0}=\sqrt{S_{s}(\mathbf{q}=0)/N}$ [see Eq. (34)]. For infinitesimal $H$, the magnetization becomes nonzero and continuously increases with increasing $H$ irrespective of $K$. This indicates continuous phase transitions at $H=0$. With further increasing $H$, the magnetization shows a continuous change between the double-$Q^{\prime}$ chiral stripe I and the single-$Q$ conical state for small $K$ and between the double-$Q^{\prime}$ chiral stripes I, II and the double-$Q$ state for large $K$. The result indicates that the transition is of second-order, although a more careful finite-size scaling is required to settle this point, especially, for large $K$. While further increasing $H$, the magnetization in both single-$Q$ conical and double-$Q$ states smoothly connect to the saturation value $M=1$ in the fully-polarized state at $H=2$.

VI.2 Triangular lattice

Figure 16 shows the $H$-$K$ phase diagram for the triangular lattice case, which consists of six phases in addition to the fully-polarized state for $H\geq 2$. We describe the details of each phase in what follows. The real-space spin configuration and the spin and chirality structure factors for each phase are shown in Fig. 17.

triple-$Q^{\prime}$ chiral stripe I (3$Q^{\prime}$-I)

This phase occupies a low-$H$ and low-$K$ region, appearing next to the double-$Q^{\prime}$ chiral stripe at $H=0$. In this phase, the $z$ component of the spin structure factor has a single-$Q$ peak at $\mathbf{Q}_{3}$, while the $xy$ component exhibits two peaks at $\mathbf{Q}_{1}$ and $\mathbf{Q}_{2}$ with different intensities, as shown in Fig. 17(a). Meanwhile, the chirality structure factor shows a single-$Q$ peak at $\mathbf{Q}_{3}$.

triple-$Q$ $n_{\rm sk}=2$ skyrmion crystal (3$Q$ with $n_{\rm sk}=2$)

This phase occupies a higher-$K$ region next to the triple-$Q^{\prime}$ chiral stripe I state at low $H$, evolving from the $H=0$ state with spin canting in the field direction. The $xy$ component of the spin structure factor exhibits two dominant peaks and a single subdominant peak, while the $z$ component shows a single dominant peak and two subdominant peaks, as shown in Fig. 17(b). The amplitudes of the total spin structure factor at $\mathbf{Q}_{1}$, $\mathbf{Q}_{2}$, and $\mathbf{Q}_{3}$ are equivalent to each other.

triple-$Q^{\prime}$ chiral stripe II (3$Q^{\prime}$-II)

This phase occupies a slightly high-$H$ region of the triple-$Q$ $n_{\rm sk}=2$ skyrmion crystal, as shown in Fig. 16. The spin configuration shown in Fig. 17(c) resembles that in the triple-$Q^{\prime}$ chiral stripe I phase shown in Fig. 17(a): it is characterized by the double-$Q$ modulation in the $xy$ component and the single-$Q$ modulation in the $z$ component of the spin structure factor. The difference is in the peak intensities of the $xy$ component: the two peaks have equal weights in the triple-$Q^{\prime}$ chiral stripe II phase, while they are different in the triple-$Q^{\prime}$ chiral stripe I phase. In the entire region of this phase, the amplitudes of the two peaks in the $xy$ component are almost the same as that of the single peak in the $z$ component, but the latter is slightly larger than the former.

triple-$Q$ $n_{\rm sk}=1$ skyrmion crystal (3$Q$ with $n_{\rm sk}=1$)

This phase is another type of crystallization of skyrmions, which has been found in chiral and frustrated magnets Mühlbauer et al. (2009); Yu et al. (2010); Nagaosa and Tokura (2013); Okubo et al. (2012); Leonov and Mostovoy (2015); Lin and Hayami (2016); Hayami et al. (2016a). It occupies a large portion of the intermediate-$H$ region, for a slightly larger $H$ than the triple-$Q^{\prime}$ chiral stripe II. Figure 17(d) shows a typical real-space spin configuration in the skyrmion crystal. The skyrmion cores drawn by the blue regions form a triangular lattice with the lattice constant $4\pi/(\sqrt{3}|\mathbf{Q}_{\nu}|)$. The spin structure factors have six peaks with equal intensities, as shown in Fig. 17(d), which indicates the presence of the $C_{6}$ symmetry. This phase also accompanies a uniform net scalar chirality, which is around half of that in the triple-$Q$ $n_{\rm sk}=2$ skyrmion crystal; this phase is characterized by the topological number $n_{\rm sk}=1$. The remarkable point of this skyrmion crystal is that it is stable even at $T=0$. This is in contrast to the similar state obtained for a localized spin model with isotropic Heisenberg interactions Okubo et al. (2012), which exists only at nonzero temperature and is replaced with the single-$Q$ conical ordering at low temperature. It is also worth noting that the skyrmion crystal state obtained here shows the degeneracy with respect to the chiral symmetry as well as the rotational symmetry along the $z$ axis, in contrast to the absence in the Heisenberg model with the DM interaction Lin and Hayami (2016).

triple-$Q^{\prime}$ chiral stripe III (3$Q^{\prime}$-III)

This phase occupies the largest portion of the phase diagram in Fig. 16, appearing next to the triple-$Q$ $n_{\rm sk}=1$ skyrmion crystal upon increasing $H$. The spin configuration resembles that in the triple-$Q^{\prime}$ chiral stripe II state, while the $xy$ component of the spin structure factor shows a small peak at $\mathbf{Q}_{1}-\mathbf{Q}_{2}$, as shown in Fig. 17(e). This phase is similar to that found in the Kondo lattice model at high field in Ref. Ozawa et al., 2016b.

triple-$Q^{\prime}$ chiral stripe IV (3$Q^{\prime}$-IV)

This phase occupies a high-$H$ and high-$K$ region of the phase diagram in Fig. 16, next to the triple-$Q^{\prime}$ chiral stripe III. In this phase, the spin configuration is also characterized by the triple-$Q$ modulation: the $xy$ component of the spin structure factor has two dominant peaks with additional four subdominant peaks, while the $z$ component shows four dominant peaks, as shown in Fig. 17(f).

Magnetization curve and net scalar chirality

Figure 18(a) shows the magnetization curve at several $K$. As in the square lattice case in Fig. 15, the magnetization becomes nonzero for infinitesimal $H$ and continuously increases with increasing $H$ irrespective of $K$, which indicates continuous phase transition at $H=0$. With further increasing $H$, the magnetization jump appears between the triple-$Q$ $n_{\rm sk}=2$ skyrmion crystal and the triple-$Q^{\prime}$ chiral stripe II. The similar magnetization jumps are found in phase boundaries between the triple-$Q$ $n_{\rm sk}=1$ skyrmion crystal and the triple-$Q^{\prime}$ chiral stripe II, III. Thus, for $K=0.1$, two jumps in the magnetization are found at $H\simeq 0.5$ and $1.0$, whereas three jumps are found for $K=0.4$ and 0.8, at $H\simeq 0.15$, $0.30$, $0.95$ and $H\simeq 0.15$, $0.35$, $0.95$, respectively, as shown in Fig. 18(a). Other phase transitions, e.g., between triple-$Q^{\prime}$ chiral stripe and fully-polarized state are continuous.

On the other hand, Fig. 18(b) shows $H$ dependence of the net scalar chirality at several $K$. The net scalar chirality is given by $\chi^{\rm sc}=(1/N)\sum_{\mathbf{R}\in\mu=(u,d)}\chi^{\rm sc}_{\mathbf{R}}$, where $\mu=(u,d)$ represent upward and downward triangles, respectively. As shown in Fig. 18(b), there are two phases which show the net scalar chirality, the triple-$Q$ $n_{sk}=2$ skyrmion crystal stabilized in the low $H$ region and the triple-$Q$ $n_{sk}=1$ skyrmion crystal in the intermediate $H$ region. In both phases, the scalar chirality is almost independent of $H$, while it changes discontinuously at their phase boundaries.

VI.3 Comparison to the Kondo lattice model

Let us compare our results with the previous studies in the Kondo lattice model. In Ref. Ozawa et al., 2016b, the Kondo lattice model on the triangular lattice was studied in an applied magnetic field. There found the phase sequence from the triple-$Q$ $n_{\rm sk}=2$ skyrmion crystal, the triple-$Q$ $n_{\rm sk}=1$ skyrmion crystal, the triple-$Q^{\prime}$ chiral stripe III (3$Q^{\prime}$-III), and to the fully-polarized state with increasing $H$. Similar sequence is found in our result in the region $0.37\lesssim K\lesssim 0.53$ in Fig. 16. The only difference is that the triple-$Q$ chiral stripe II appears between the triple-$Q^{\prime}$ $n_{\rm sk}=2$ skyrmion crystal and the triple-$Q$ $n_{\rm sk}=1$ skyrmion crystal in our results. A possible reason for this discrepancy is due to the lack of the effect of the magnetic field on Green’s function in our effective model, which gives rise to the spin dependent interaction via the Zeeman coupling in Eq. (65). Among them, the third-order contributions with respect to $J$ arise in nonzero $H$, and may prefer the triple-$Q$ skyrmion crystals to the triple-$Q$ chiral stripe; see Appendix A. Another possibility is the charge degree of freedom in the Kondo lattice model. In our effective spin model, it is implicitly assumed that the electron density is unchanged while changing the magnetic field, while the results in the Kondo lattice model were calculated for a fixed chemical potential Ozawa et al. (2016b). Such field dependence of the electron density may result in the discrepancy.

Meanwhile, for the square lattice case, the authors have found the phase transition from the double-$Q^{\prime}$ chiral stripe to the double-$Q$ state in the Kondo lattice model with increasing $H$ in Ref. Oza, . This corresponds to the region for $0.27\lesssim K$ in our effective spin model. This agreement supports the validity of our effective model.

Although the detailed comparison is left for future study, we stress that our effective model reproduces well the exotic magnetic phases found in the Kondo lattice model. The effective spin model will be useful to explore further exotic phases by smaller computational efforts compared to those for itinerant electron models.

VII Summary and Concluding Remarks

To summarize, we have constructed an effective spin model for describing magnetic instabilities in itinerant magnets. Taking one of the fundamental models for itinerant magnets, the Kondo lattice model, and carefully examining the spin scattering processes by the perturbation in terms of the spin-charge coupling, we proposed a simple model including the bilinear and biquadratic interactions with particular wave numbers dictated by the Fermi surface. The bilinear interaction has the same form as the RKKY interaction, and prefers a helical magnetic order specified by a single wave number. On the other hand, the biquadratic interaction, which is deduced from the dominant contributions in the higher-order perturbations, causes an instability toward noncollinear or noncoplanar ordering specified by multiple wave numbers. We have tested the validity of the effective model by calculating the ground-state phase diagrams by Monte Carlo simulation and comparing the results with those for the Kondo lattice model. The comparison on square and triangular lattices shows a good agreement with the previous findings in the Kondo lattice model Ozawa et al. (2016a, b): we have demonstrated that our effective spin model reproduces the noncoplanar double-$Q$ instability in both square and triangular lattice cases and the triple-$Q$ skyrmion with topological number of two in the triangular lattice case. The good agreement indicates that our effective model captures the underlying mechanism of the magnetic instabilities toward noncollinear and noncoplanar orderings in itinerant magnets: the key ingredient is the competition between bilinear and biquadratic interactions in momentum space.

We have also extended our study by introducing the external magnetic field to our effective spin model. The simplicity of the model allows us to systematically investigate the magnetic phase diagram by much smaller computational costs compared to the original Kondo lattice model. In addition to the phases already found in the previous studies for the Kondo lattice model Ozawa et al. (2016b), we found several new phases in an applied magnetic field. The results demonstrate the efficiency of our model for further exploration of exotic magnetic phases in itinerant magnets.

Our effective spin model is constructed on the basis of a fundamental property of the spin-charge system: multiple peaks in the bare susceptibility whose wave numbers are related with the lattice symmetry. As this is commonly seen in all the centrosymmetric lattices in two and three dimensions, we believe that our model is applicable to a wide class of itinerant magnets. In particular, the instability toward noncollinear and noncoplanar orderings induced by the effective biquadratic interaction will be a universal feature irrespective of the details of the system. Also, we note that our effective spin model can be easily extended to more complicated situations, by including, e.g., anisotropic interactions, single-ion anisotropy, and dipole-dipole interactions. Such extensions will be helpful to investigate and revisit unconventional magnetic phase diagrams from the viewpoint beyond the RKKY mechanism.

Finally, let us comment on candidate materials, whose magnetic properties might be described by our spin model. One of the candidates is a scandium thiospinel MnSc2S4, which has been recently identified to show a triple-$Q$ vortex crystal under an external magnetic field Gao et al. (2016). Another candidate is the strontium iron perovskite oxide SrFeO3, which shows several phases depending on magnetic field and temperature Ishiwata et al. (2011). Interestingly, the transport measurements imply that the phases include multiple-$Q$ states in spite of the negligibly small contribution of the spin-orbit coupling. Our model will make a significant contribution to understand these exotic magnetisms. The detailed comparison will be left for future study.

Acknowledgements.

R.O. is supported by the Japan Society for the Promotion of Science through a research fellowship for young scientists and the Program for Leading Graduate Schools (ALPS). This research was supported by KAKENHI (No. 16H06590), the Strategic Programs for Innovative Research (SPIRE), MEXT, and the Computational Materials Science Initiative (CMSI), Japan.

Appendix A Third-order effective magnetic interactions in a magnetic field

In this appendix, we show the expression of the third-order magnetic interactions, which appears under the external magnetic field. By substituting $\mathbf{q}=\mathbf{Q}_{1}$, $\mathbf{Q}_{2}$, and $\mathbf{Q}_{3}$ in the expansion of the free energy in Eq. (5), we obtain the third-order contributions, which are given by

[TABLE]

where H.c. represent the Hermitian conjugate. The coefficients are represented by

[TABLE]

where $C_{1}$ to $C_{6}$ represent the scattering processes by a single-$Q$ wave number, while $D_{1}$ to $D_{4}$ the scattering processes by multiple-$Q$ wave numbers. We may need to take into account these terms in the effective spin model in nonzero field, which is left for future study.