Magnetic hedgehog lattices in noncentrosymmetric metals

TL;DR

This paper investigates the microscopic origins and stability of magnetic hedgehog lattices in noncentrosymmetric metals, revealing their stabilization mechanisms and phase transitions under magnetic fields through theoretical modeling.

Contribution

It demonstrates that HLs are stabilized at zero magnetic field by spin-orbit and spin-charge interactions, providing a microscopic understanding of their origin and behavior.

Findings

HLs are stabilized at zero field by spin interactions.

Phase diagram includes topological transitions with monopole pair annihilation.

Different types of HLs can exist even at zero magnetic field.

Abstract

The magnetic hedgehog lattice (HL) is a noncoplanar magnetic texture with a periodic array of magnetic monopoles and anti-monopoles. Despite phenomenological and numerical studies thus far, there remain open issues on the microscopic origin, especially with respect to the recent experimental findings of two different types of HLs even at zero magnetic field. Here, we study the stability of the HLs for an effective spin model with long-range interactions arising from itinerant nature of electrons. By variational calculations and simulated annealing, we find that the HLs are stabilized in the ground state at zero magnetic field by the synergetic effect of the anti-symmetric exchange interactions generated by the spin-orbit coupling and the multiple-spin interactions generated by the spin-charge coupling. We also clarify the phase diagram in the magnetic fields, which includes topological phase transitions with pair annihilation of the monopoles and anti-monopoles depending on the field directions.

Equations

$$\mathcal{H}=$$

$$+\sum _{\mathbf{k}\sigma {\sigma }^{\prime }}{\mathbf{g}}_{\mathbf{k}}\cdot c_{\mathbf{k}\sigma }^{†}{\mathbf{\sigma }}_{\sigma {\sigma }^{\prime }}{c}_{\mathbf{k}{\sigma }^{\prime }}^{},$$

$$\displaystyle\mathcal{H}=\sum_{\eta}\Big{[}$$

$$\displaystyle-i{\mathbf{D}}_{\eta}\cdot{\mathbf{S}}_{\mathbf{Q}_{\eta}}\times{\mathbf{S}}_{-\mathbf{Q}_{\eta}}\Big{]}-\sum_{l}\mathbf{h}\cdot\mathbf{S}_{\mathbf{r}_{l}}.$$

$${\mathbf{S}}_{{\mathbf{r}}_l}\propto \sum _{\eta =1}^n{\hat{\mathbf{a}}}_{\eta }\cos {\mathcal{Q}}_{\eta l},$$

$${\mathbf{S}}_{{\mathbf{r}}_l}\propto \sum _{\eta =1}^n({\hat{\mathbf{b}}}_{\eta }\sin {\mathcal{Q}}_{\eta l}+{\hat{\mathbf{c}}}_{\eta }\cos {\mathcal{Q}}_{\eta l}),$$

$${\mathbf{S}}_{{\mathbf{r}}_l}\propto \sqrt{1-u^2}{\hat{\mathbf{b}}}_1\sin {\mathcal{Q}}_{1l}+\sqrt{1-u^2}{\hat{\mathbf{c}}}_1\cos {\mathcal{Q}}_{1l}+u{\hat{\mathbf{a}}}_1,$$

$$m=\frac{1}{N}\sum _l{\mathbf{S}}_{{\mathbf{r}}_l}\cdot \hat{\mathbf{h}}.$$

$$\chi =\frac{m(h+\Delta h)-m(h)}{\Delta h}.$$

$$m_{\mathbf{q}}=\sqrt{\frac{S(\mathbf{q})}{N}},$$

$$S(\mathbf{q})=\frac{1}{N}\sum _{l,l^{\prime }}{\mathbf{S}}_{{\mathbf{r}}_l}\cdot {\mathbf{S}}_{{\mathbf{r}}_{l^{\prime }}}{e}^{i\mathbf{q}\cdot ({\mathbf{r}}_l-{\mathbf{r}}_{l^{\prime }})}.$$

$$Q_{\mathrm{m}}({\mathbf{r}}_{\mathrm{c}})=\frac{1}{4\pi }\sum _{p\in \mathrm{unit}\mathrm{cube}}{\mathbf{\Omega }}_p\cdot {\hat{\mathbf{n}}}_p,$$

$${\Omega }_i=2{\tan }^{-1}\left\{\frac{{\mathbf{S}}_1\cdot \left({\mathbf{S}}_2\times {\mathbf{S}}_3\right)}{1+{\mathbf{S}}_1\cdot {\mathbf{S}}_2+{\mathbf{S}}_2\cdot {\mathbf{S}}_3+{\mathbf{S}}_3\cdot {\mathbf{S}}_1}\right\},$$

$${\mathbf{\Omega }}_p=\sum _{i\in p}{\Omega }_i{\hat{\mathbf{n}}}_p.$$

$$N_{\mathrm{m}}=\sum _{{\mathbf{r}}_{\mathrm{c}}\in V_{\mathrm{m}}}|Q_{\mathrm{m}}({\mathbf{r}}_{\mathrm{c}})|,$$

$$d_{\mathrm{m}}=\min |{\mathbf{r}}_{\mathrm{c}}^{\mathrm{m}}-{\mathbf{r}}_{\mathrm{c}}^{\mathrm{a}}|,$$

$${\chi }_{\mathrm{sc}}^{\gamma }({\mathbf{r}}_l)=\frac{1}{2}\sum _{\alpha \beta {\nu }_{\alpha }{\nu }_{\beta }}{ϵ}^{\alpha \beta \gamma }{\nu }_{\alpha }{\nu }_{\beta }{\mathbf{S}}_{{\mathbf{r}}_l}\cdot ({\mathbf{S}}_{{\mathbf{r}}_l+{\nu }_{\alpha }{\hat{\mathbf{\delta }}}_{\alpha }}\times {\mathbf{S}}_{{\mathbf{r}}_l+{\nu }_{\beta }{\hat{\mathbf{\delta }}}_{\beta }}),$$

$${\chi }_{\mathrm{sc}}=\frac{1}{N}\sum _{\gamma l}{\chi }_{\mathrm{sc}}^{\gamma }({\mathbf{r}}_l).$$

$$S^{\parallel }(\mathbf{q})=\frac{1}{2}\{S^{xx}(\mathbf{q})+S^{yy}(\mathbf{q})\}+S^{xy}(\mathbf{q}),$$

$$S^{\mu \nu }(\mathbf{q})=\frac{1}{N}\sum _{l,l^{\prime }}{S}_{{\mathbf{r}}_l}^{\mu }{S}_{{\mathbf{r}}_{l^{\prime }}}^{\nu }{e}^{i\mathbf{q}\cdot ({\mathbf{r}}_l-{\mathbf{r}}_{l^{\prime }})}.$$

$${\chi }_{\mathrm{sc}}^{\perp }(\mathbf{q})=\frac{1}{2}\{{\chi }_{\mathrm{sc}}^{xx}(\mathbf{q})+{\chi }_{\mathrm{sc}}^{yy}(\mathbf{q})\}+{\chi }_{\mathrm{sc}}^{xy}(\mathbf{q}),$$

$${\chi }_{\mathrm{sc}}^{\mu \nu }(\mathbf{q})=\frac{1}{N}\sum _{l,l^{\prime }}{\chi }_{\mathrm{sc}}^{\mu }({\mathbf{r}}_l){\chi }_{\mathrm{sc}}^{\nu }({\mathbf{r}}_{l^{\prime }})e^{i\mathbf{q}\cdot ({\mathbf{r}}_l-{\mathbf{r}}_{l^{\prime }})}.$$

$$S^{\perp }(\mathbf{q})=\frac{2}{3}\{S(\mathbf{q})-S^{xy}(\mathbf{q})-S^{yz}(\mathbf{q})-S^{zx}(\mathbf{q})\}.$$

Full text

I Introduction

Chirality, often termed as handedness, is a key concept in a broad field of science, ranging from particle physics to biology. In condensed matter physics, chiral magnetic textures, which break both inversion and mirror symmetries in addition to time-reversal symmetry, have recently attracted considerable attention for potential applications to next-generation electronic devices. There are a variety of the chiral magnetic textures, such as skyrmion lattices Nagaosa and Tokura (2013) and chiral soliton lattices Togawa et al. (2016). Noncollinear and noncoplanar spin arrangements in these textures generate emergent electromagnetic fields through the Berry phase mechanism, which induce unconventional transport, optical, and magnetoelectric properties Tokura and Seki (2010); Mochizuki and Seki (2015); Tokura and Nagaosa (2018).

Recently, a three-dimensional chiral magnetic texture, which is called the hedgehog lattice (HL), was discovered in the $B$20-type compound MnGe Tanigaki et al. (2015); Kanazawa et al. (2017). The magnetic structure is characterized by cubic three wave vectors, and hence, it is referred as the triple-$Q$ hedgehog lattice ($3Q$-HL) [Fig. 1(b)]. The $3Q$-HL has a periodic array of hyperbolic hedgehog and anti-hedgehog spin textures, which generates an emergent magnetic field with a periodic array of radial hedgehogs and anti-hedgehogs regarded as magnetic monopoles and anti-monopoles, as shown in Fig. 1(c) Kanazawa et al. (2012, 2016); Zhang et al. (2016). The peculiar magnetic field was discussed as a source of the enormous topological Hall effect Kanazawa et al. (2011) and thermoelectric effect Shiomi et al. (2013); Fujishiro et al. (2018). In addition, by a substitution of Ge by Si, the $3Q$-HL changes into a different HL characterized by tetrahedral four wave vectors, dubbed the quadruple-$Q$ hedgehog lattice ($4Q$-HL) [Fig. 1(a)] Fujishiro et al. (2019). Remarkably, the magnetic periods of these $3Q$- and $4Q$-HLs are very short $\sim 2$-$3$ nm, in contrast to most of the skyrmion lattices.

Such magnetic HLs have been theoretically studied prior to the experimental discovery, e.g., by the Ginzburg-Landau theory Binz and Vishwanath (2006), variational calculations Park and Han (2011), and Monte Carlo (MC) simulations Yang et al. (2016). The variational study for a classical spin model showed that the $3Q$-HL is not stabilized, whereas the $4Q$-HL is obtained in an applied magnetic field Park and Han (2011). The $4Q$-HL in a field was also confirmed by MC simulations Yang et al. (2016). The previous studies, however, do not predict the stable HLs in the absence of magnetic fields, contradicting the experimental observations. Furthermore, to account for the short-period twist, the localized spin picture requires a large Dzyaloshinskii-Moriya (DM) interaction Dzyaloshinsky (1958); Moriya (1960), but it was estimated to be very weak Gayles et al. (2015); Koretsune et al. (2015); Kikuchi et al. (2016). Indeed, recent analyses based on first-principles calculations showed that the stable HLs are not obtained by two spin interactions including the DM interaction Grytsiuk et al. (2019). The importance of four- and six-spin interactions including spin chirality was also proposed Brinker et al. (2019); Grytsiuk et al. (2020).

In this paper, we study the stability of $4Q$- and $3Q$-HLs from a different viewpoint from the previous studies, by taking into account itinerant nature of electrons. We consider an effective model with long-range exchange interactions originating from the coupling between charge, spin, and orbital degrees of freedom. By variational calculations and simulated annealing, we show that the model realizes both $4Q$- and $3Q$-HLs at zero field, through the cooperation between the DM-type asymmetric exchange interactions arising from the spin-orbit coupling and the multiple-spin interactions from the spin-charge coupling. We also study the effect of an applied magnetic field on these HLs. Depending on the field directions, we find that the system exhibits multiple phase transitions while changing from the $4Q$- and $3Q$-HLs to the forced ferromagnetic (FFM) state. Notably, we show that some of them are topological phase transitions with pair annihilation of the monopoles and anti-monopoles. We demonstrate how the pair annihilation takes place by tracing the positions of the monopoles and anti-monopoles.

The rest of the paper is organized as follows. In Sec. II, we introduce the effective spin model derived from an itinerant electron model. In Sec. III, we describe the methods that we use in this study to investigate the ground state of the effective spin model. In Sec. IV, we show the phase diagram at zero field including the HLs. In Sec. V, we show the phase diagram in magnetic fields applied in three symmetric directions. In Sec. VI, we discuss field-induced topological phase transitions caused by pair annihilation of monopoles and anti-monopoles. Section VII is devoted to the summary.

II Model

In this section, we present the model which we use in the present study. Starting from an itinerant electron model with spin-charge and spin orbit couplings in Sec. II.1, we discuss the effective model with long-ranged exchange interactions induced by the itinerant nature of electrons in Sec. II.2.

II.1 Itinerant electron model

In order to investigate the microscopic origin of magnetic HLs, we begin with a minimal model including itinerant electrons, an extended Kondo lattice model that describes the coupling between the itinerant electron spins and localized magnetic moments. While the Kondo lattice model has been studied for $f$ electron systems, where the $f$ electrons comprises the localized moments Stewart (1984); Gegenwart et al. (2008), we note that it is also regarded as an effective model for the Hubbard-type models, which have been used widely, e.g., for $d$ electron systems, at the level of the mean-field approximation Martin and Batista (2008). In the current study, we include an anti-symmetric spin-orbit coupling arising from spatial inversion symmetry breaking in noncentrosymmetric systems. The Hamiltonian in the wave-number representation is given by

[TABLE]

where $c^{\dagger}_{\mathbf{k}\sigma}$ ($c^{\;}_{\mathbf{k}\sigma}$) is a creation (annihilation) operator of an itinerant electron with wave vector $\mathbf{k}$ and spin $\sigma=\uparrow$ or $\downarrow$. The first term describes the kinetic energy of itinerant electrons; $\varepsilon_{\mathbf{k}}$ is the energy dispersion and $\mu$ is the chemical potential. The second term is for the Kondo coupling between itinerant electron spins and localized spin moments; ${\bm{\sigma}}=(\sigma^{x},\sigma^{y},\sigma^{z})$ is the vector of Pauli matrices, and $\mathbf{S}_{\mathbf{q}}=\frac{1}{\sqrt{N}}\sum_{l}\mathbf{S}_{\mathbf{r}_{l}}e^{-i\mathbf{q}\cdot\mathbf{r}_{l}}$ is the Fourier transform of a localized moment $\mathbf{S}_{\mathbf{r}_{l}}=(S^{x}_{\mathbf{r}_{l}},S^{y}_{\mathbf{r}_{l}},S^{z}_{\mathbf{r}_{l}})$ defined at site $l$, where $N$ is the number of lattice sites. For simplicity, $\mathbf{S}_{\mathbf{r}_{l}}$ is regarded as a classical spin with the length $|\mathbf{S}_{\mathbf{r}_{l}}|=1$. $J_{\mathrm{K}}$ is the exchange coupling constant whose sign is irrelevant for the classical spins. The last term represents the anti-symmetric spin-orbit coupling induced by spatial inversion symmetry breaking; $\mathbf{g}_{\mathbf{k}}=(g^{x}_{\mathbf{k}},g^{y}_{\mathbf{k}},g^{z}_{\mathbf{k}})$ is called the g-vector, which plays an important role in chiral magnets. In the following, we consider the model on a simple cubic lattice with the lattice constant being unity for simplicity; noncentrosymmetric nature is effectively taken into account in the g-vector $\mathbf{g}_{\mathbf{k}}$ with an odd-function of $\mathbf{k}$.

II.2 Effective spin model

In general, the coupling between itinerant electrons and localized spins generates effective exchange interactions between the localized spins. For instance, in the strong-coupling case with $J_{\mathrm{K}}\gg|\varepsilon_{\mathbf{k}}-\mu|$, an effective ferromagnetic interaction is generated to maximize the kinetic energy of itinerant electrons by aligning neighboring spins, which is called the double-exchange interaction Zener (1951); Anderson and Hasegawa (1955). On the other hand, in the weak-coupling case with $J_{\mathrm{K}}\ll|\varepsilon_{\mathbf{k}}-\mu|$, the effective magnetic interaction becomes long-ranged and oscillating in space, which is called the Ruderman-Kittel-Kasuya-Yosida (RKKY) interaction Ruderman and Kittel (1954); Kasuya (1956); Yosida (1957). In this study, we consider the weak-coupling case of the model in Eq. (1) by an effective spin model derived by perturbation expansion in terms of $J_{\mathrm{K}}$. Our model includes a higher-order effect of the spin-charge coupling beyond the RKKY interaction discussed in the previous studies Akagi et al. (2012); Hayami and Motome (2014); Hayami et al. (2017), and also a DM-type interaction originating from the spin-orbit coupling in the last term in Eq. (1) Hayami and Motome (2018). The Hamiltonian reads

[TABLE]

The first term denotes the RKKY interaction, which is derived by the second-order perturbation with respect to $J_{\mathrm{K}}$ Ruderman and Kittel (1954); Kasuya (1956); Yosida (1957). In general, this tends to stabilize a spiral magnetic texture. The second term is the biquadratic interaction, which is most relevant among the higher-order perturbations with respect to $J_{\mathrm{K}}$ Hayami et al. (2017). Hereafter, we consider the positive coupling constant $K>0$, which is known to prefer noncollinear and noncoplanar spin configurations Akagi et al. (2012); Hayami and Motome (2014); Hayami et al. (2017). The third term represents a DM-type interaction arising from the anti-symmetric spin-orbit coupling, which is derived by the second-order perturbation with respect to $J_{\mathrm{K}}$ Hayami and Motome (2018). This also brings a twist in spin textures, and plays a role in not only choosing the chirality but also giving an anisotropy in spin space. Note that we ignore other anisotropic exchange interactions originating from the anti-symmetric spin-orbit coupling, for simplicity Hayami and Motome (2018). The last term describes the Zeeman coupling to an external magnetic field $\mathbf{h}$.

In Eq. (2), all the exchange interactions are long-ranged in real space and specified by particular wave numbers $\mathbf{Q}_{\eta}$. This inherits the itinerant nature of electrons; specifically, the wave vectors $\mathbf{Q}_{\eta}$ are set by the multiple maxima in the spin-dependent bare susceptibility of itinerant electrons Hayami and Motome (2014); Hayami et al. (2017). Corresponding to the $3Q$- and $4Q$-HLs, we assume two sets of $\mathbf{Q}_{\eta}$: One is a set of the tetrahedral wave vectors as $\mathbf{Q}_{1}=(Q,-Q,-Q)$, $\mathbf{Q}_{2}=(-Q,Q,-Q)$, $\mathbf{Q}_{3}=(-Q,-Q,Q)$, and $\mathbf{Q}_{4}=(Q,Q,Q)$ [Fig. 1(a)], and the other is a set of the cubic wave vectors as $\mathbf{Q}_{1}=(Q,0,0)$, $\mathbf{Q}_{2}=(0,Q,0)$, and $\mathbf{Q}_{3}=(0,0,Q)$, which are orthogonal to each other [Fig. 1(b)]. In the following calculations, we set $Q=\pi/4$ (period of eight lattice sites); we confirm that the following results remain qualitatively the same for different choices of $Q$. Although the direction of $\mathbf{D}_{\eta}$ is independent of that of $\mathbf{Q}_{\eta}$ in general, we assume $\mathbf{D}_{\eta}\parallel\mathbf{Q}_{\eta}$ that stabilizes proper-screw type spin textures 111The perturbation expansion for Eq. (1) leads to $\mathbf{D}_{\eta}\parallel\mathbf{Q}_{\eta}$ Hayami and Motome (2018). We note that the HLs can be composed of superpositions of the proper screws. The magnetic field $\mathbf{h}$ is applied along the [001], [110], and [111] directions as shown in the right panels of Figs. 1(a) and 1(b). We set the energy scale as $J=1$. We consider the system with $N=16^{3}$ spins under periodic boundary conditions. We confirmed that the following results remain the same for $N=24^{3}$ spins (not shown here).

III Method

In this section, we present the methods to study the ground state of the model in Eq. (2). At zero magnetic field, we mainly adopt variational calculations by comparing the energy of several different spin states, as introduced in Sec. III.1. In addition, we use simulated annealing, which is introduced in Sec. III.2, not only to confirm the variational results but also to study the ground state in an applied magnetic field where it is difficult to infer the variational states.

III.1 Variational calculations

In the variational calculations, we consider the following spin textures as the variational states at zero magnetic field. The simplest one is given by

[TABLE]

where $\hat{\mathbf{a}}_{\eta}$ is the unit vector parallel to $\mathbf{Q}_{\eta}$ and $\mathcal{Q}_{\eta l}=\mathbf{Q}_{\eta}\cdot\mathbf{r}_{l}+\varphi_{\eta}$ ($\varphi_{\eta}$ represents the phase shift); $n=1,2,3$ for the $3Q$ case and $n=1,2,3,4$ for the $4Q$ case. This is a set of nonchiral states that has no energy gain from the DM-type interaction. Another variational state is a chiral one described as the equal superpositions of proper screws,

[TABLE]

where $\hat{\mathbf{b}}_{\eta}$ and $\hat{\mathbf{c}}_{\eta}$ are the unit vectors orthogonal to $\hat{\mathbf{a}}_{\eta}$ and each other ($\hat{\mathbf{a}}_{\eta}$, $\hat{\mathbf{b}}_{\eta}$, and $\hat{\mathbf{c}}_{\eta}$ form a right-handed system). Note that the $n=3$ ($n=4$) state for $3Q$ ($4Q$) corresponds to the $3Q$($4Q$)-HL shown in Fig. 1(a)[(b)]. In addition, we include another variational state called the double-$Q$ chiral stripe (2$Q$-CS) found in the previous study Ozawa et al. (2016),

[TABLE]

where $u=v\sin\mathcal{Q}_{2l}$. In the variational calculations, we compare the energy for all the variational states by varying $\varphi_{\eta}$ from [math] to $Q$ and $v$ from [math] to $1$ to find the lowest-energy candidate for the ground state.

III.2 Simulated annealing

In the simulated annealing, we numerically find the candidate for the ground state by mean of MC simulation. We gradually reduce the temperature of the system from $T=1$ to $T=10^{-5}$ with a condition $T_{n}=10^{-0.1n}$, where $T_{n}$ is the temperature in the $n$th step. During the annealing, we spend a total of $10^{5}-10^{6}$ MC sweeps by using the standard Metropolis algorithm. After annealing at a particular value of the field strength $h=|\mathbf{h}|$, we increase or decrease $h$ successively by $\Delta h=0.01$. At every shift by $\Delta h$, we heat the system up to $T=10^{-3}$ and cool down again to $T=10^{-5}$ by annealing. Carefully comparing the energy by starting from various values of $h$, we map out the magnetic phase diagram.

For the state obtained by the simulated annealing, we calculate the magnetization per site along the field direction,

[TABLE]

where $\hat{\mathbf{h}}$ is the unit vector in the field direction, and the magnetic susceptibility,

[TABLE]

To identify the multiple-$Q$ magnetic orders, we also calculate the magnetic moment with wave vector $\mathbf{q}$,

[TABLE]

where $S(\mathbf{q})$ is the spin structure factor defined by

[TABLE]

In addition, following Ref. Yang et al. (2016), we define the monopole charge in each unit cube by using the fluxes $\bm{\Omega}_{p}$ penetrating six square plaquettes of the cube as Okumura et al.

[TABLE]

where $\mathbf{r}_{\mathrm{c}}$ is the center position of the unit cube and $\hat{\mathbf{n}}_{p}$ is the normal unit vector of the $p$th plaquette pointing outward of the cube. We compute the flux $\bm{\Omega}_{p}$ by dividing the $p$th plaquette into two triangles and taking the sum of the solid angles of three spins on the two triangles $i=1$ and $2$. Each solid angle is calculated by

[TABLE]

where $\mathbf{S}_{1}$, $\mathbf{S}_{2}$, and $\mathbf{S}_{3}$ are the three spins on the $i$th triangle in the clockwise order viewed from the center of the cube, and the sign of $\Omega_{i}$ is taken to be the same as that of $\mathbf{S}_{1}\cdot\left[\mathbf{S}_{2}\times\mathbf{S}_{3}\right]$: $\Omega_{i}\in[-2\pi,2\pi]$. The flux $\bm{\Omega}_{p}$ is defined as a perpendicular vector to the $p$th plaquette as

[TABLE]

By substituting Eq. (12) into Eq. (10), we obtain the monopole charge $Q_{\mathrm{m}}(\mathbf{r}_{\mathrm{c}})$. This quantity detects the monopoles and anti-monopoles as it takes the value of +1 (-1) when a monopole (anti-monopole) exists in the unit cube. The monopoles and anti-monopoles are connected by a flow of the flux $\bm{\Omega}_{p}$ in Eq. (12). We compute the total number of monopoles and anti-monopoles in the magnetic unit cell, $N_{\mathrm{m}}$, as

[TABLE]

where $V_{\mathrm{m}}$ is the magnetic unit cell ($8^{3}$ sites in the following calculations). We also measure the distances between the monopoles and anti-monopoles by using $\mathbf{r}_{\mathrm{c}}$ where $\mathbf{Q}_{\mathrm{m}}(\mathbf{r}_{\mathrm{c}})=\pm 1$. In particular, we compute the minimum distance between the monopole and anti-monopole by

[TABLE]

where $\mathbf{r}_{\mathrm{c}}^{\mathrm{m}}$ and $\mathbf{r}_{\mathrm{c}}^{\mathrm{a}}$ denote $\mathbf{r}_{\mathrm{c}}$ for the monopoles and anti-monopoles. This is an important quantity for not only monitoring topological phase transitions by pair annihilation between monopoles and anti-monopoles but also understanding the behavior of the net scalar spin chirality introduced below. We note, however, that $\mathbf{r}_{\mathrm{c}}^{\mathrm{m(a)}}$ gives an approximate position of the (anti-)monopole core within an accuracy of the lattice constant, and $d_{\mathrm{m}}$ changes discontinuously by definition.

Finally, we calculate the net scalar spin chirality. We define the scalar spin chirality at each lattice site ${\mathbf{r}}_{l}$ by the sum of spin triple products on four triangles on the $\alpha\beta$ plane ($\alpha,\beta=x,y,z$) as Okumura et al.

[TABLE]

where $\gamma$ is the perpendicular direction to the $\alpha\beta$ plane, $\epsilon^{\alpha\beta\gamma}$ is the Levi-Civita symbol, $\nu_{\alpha(\beta)}=\pm 1$, and $\hat{\bm{\delta}}_{\alpha(\beta)}$ is the unit translation vector in the $\alpha(\beta)$ direction. By taking the sum over all the sites and three planes, we obtain the net spin chirality:

[TABLE]

Since Eqs. (11) and (15) share the spin triple products, $\chi_{\mathrm{sc}}$ is related with the (oriented) summation of the flux $\bm{\Omega}_{p}$ in Eq. (12). As mentioned above, the flows of the flux connect the monopoles and anti-monopoles, and hence, the lengths of the flux flows, which are approximately given by the distances $|\mathbf{r}_{\mathrm{c}}^{\mathrm{m}}-\mathbf{r}_{\mathrm{c}}^{\mathrm{a}}|$, affect $\chi_{\mathrm{sc}}$. We will discuss such a relation in Sec. VI.

IV Phase diagram at zero field

First, we show the results in the absence of the magnetic field obtained by the variational calculations in Sec. III.1. Figures 2(a) and 2(b) display the magnetic phase diagrams for the 4$Q$ and 3$Q$ cases, respectively, while varying $D=|\mathbf{D}_{\eta}|$ and $K$ in Eq. (2). When $K=0$, a nonzero $D$ stabilizes the chiral $1Q$ helical state ($1Q$-H), which remains stable in the small $K$ region for $D>0$ in both $4Q$ and $3Q$ cases. On the other hand, when introducing $K$ with $D=0$, the 2$Q$-CS is stabilized in both cases, but replaced by the nonchiral 4$Q$ and 3$Q$ states in the larger $K$ region. Similar sequence of the phase transitions was found in two dimensions Ozawa et al. (2016); Hayami et al. (2017). When $D$ and $K$ are both relevant, however, we find the 4$Q$- and 3$Q$-HLs in the wide parameter range, in addition to a chiral $2Q$ state in the $3Q$ case, which is a Bloch-type vortex crystal ($2Q$-VC) Hayami and Motome (2018). We confirm the stability of these HLs also by the simulated annealing in Sec. III.2; typical spin configurations for the 4$Q$- and 3$Q$-HLs are presented in Figs. 1(a) and 1(b), respectively.

Thus, our results indicate that the $4Q$- and $3Q$-HLs are stabilized by cooperation between the RKKY interaction, the biquadratic interaction, and the DM-type interaction. In other words, both spin-charge and spin-orbit couplings play a crucial role in the stabilization of the $4Q$- and $3Q$-HLs.

From the variational calculations, we find that the stable positions of all the monopoles and anti-monopoles of the $4Q$- and $3Q$-HLs locate not at the lattice sites but at the interstitial positions (centers of unit cubes). This is concluded for the $3Q$-HL by that the optimized phase shift in Eq. (4) always takes $\varphi_{\eta}=\pi/8$. In this case, the eight spins surrounding the (anti-)monopole comprise a hyperbolic (anti-)hedgehog whose north and south poles are in the [111] direction, as shown in Fig. 1(c). Meanwhile, for the $4Q$-HL, the set of $\varphi_{\eta}$ depends on $D$ and $K$ since the four ordering vectors $\mathbb{Q}_{\eta}$ are dependent on each other. In this case, however, the eight spins comprise a hyperbolic (anti-)hedgehog with the north and south poles in the [001] direction. In both cases, the (anti-)hedgehogs generates an effective (anti-)monopole field, as shown in Fig. 1(c). We deduce that the stable monopoles and anti-monopoles centered at the interstitial positions might be ubiquitous to the systems with fixed spin length on the discrete lattice since their cores are singular points where the spins vanish in the continuum limit.

V phase transitions in magnetic fields

Next, we show the results for the phase diagrams of the model in Eq. (2) in the magnetic fields along the [001], [110], and [111] directions obtained by the simulated annealing in Sec. III.2. In Secs. V.1 and V.2, we present the results for the $4Q$ and $3Q$ cases, respectively.

V.1 4$Q$ case

Let us first discuss the $4Q$ case, whose ordering vectors are shown in Fig. 1(a). Figure 3 summarizes the results for the $4Q$-HL at $D=0.3$ and $K=0.6$.

First, we discuss the results for the [001] field, $\mathbf{h}_{001}=(0,0,h)$, shown in Fig. 3(a). As plotted in the top panel, the magnetization $m$ shows kinks at $h\simeq 0.575$, $1.395$, and $2.335$, and a small jump at $h\simeq 0.595$. Correspondingly, the magnetic susceptibility $\chi$ shows peaks at $h\simeq 0.575$ and $0.595$, a broad hump at $h\simeq 1.395$, and a shoulder at $h\simeq 2.335$. These indicate the existence of at least four phase transitions: one at $h\simeq 0.595$ is of first order, while the rest three are of second order. The magnetic moments $m_{\mathbf{Q}_{\eta}}$ plotted in the middle panel show that the four phases below $h\simeq 2.335$ are $4Q$ states with the equal amplitudes for the four $m_{\mathbf{Q}_{\eta}}$, whereas the phase for $h\gtrsim 2.335$ is a FFM state. We note that these 4$Q$ states are distinguished by the higher Fourier components of the spin structure factor $S(\mathbf{q})$ (see Appendix A).

The number of monopoles and anti-monopoles, $N_{\mathrm{m}}$, is plotted in the bottom panel. The result shows that $N_{\rm m}$ is halved and vanishes through the second-order phase transitions at $h\simeq 0.575$ and $1.395$, respectively (black dashed lines). As plotted in the inset, the minimum distance between the monopoles and anti-monopoles, $d_{\mathrm{m}}$, gets shorter from $d_{\rm m}=2\sqrt{3}$ to $1$ and $3$ to $1$ while approaching $h\simeq 0.575$ and $1.395$, respectively. These suggest that the phase transitions are topological ones caused by pair annihilation of monopoles and anti-monopoles. We will discuss the details in Sec. VI.1.

In the bottom panel, we also plot the net scalar spin chirality $\chi_{\rm sc}$, which gives rise to the topological Hall effect in itinerant electron systems Binz and Vishwanath (2008). $\chi_{\rm sc}$ rapidly increases before the phase transition at $h\simeq 0.575$. After showing a sharp peak at the phase transition at $h\simeq 0.595$, $\chi_{\rm sc}$ exhibits a broad peak at $h\sim 1$, rapidly decreases around the phase transition at $h\simeq 1.395$, and smoothly reduces to zero while approaching the phase transition to a FFM state at $h\simeq 2.335$. The change of $\chi_{\rm sc}$ is closely related with the change in the lengths of flows of the flux $\bm{\Omega}_{p}$ in Eq. (12); see Sec. VI.1.

Next, we discuss the results for the [110] field, $\mathbf{h}_{110}=\frac{1}{\sqrt{2}}(h,h,0)$, shown in Fig. 3(b). As shown in the top panel, $m$ and $\chi$ show jumps and sharp peaks, respectively, at $h\simeq 0.795$, $0.845$, and $1.435$. $m$ also has kinks at $h\simeq 1.325$ and $2.495$, where $\chi$ shows a broad hump and a shoulder, respectively. These indicate the existence of at least five phase transitions: three at $h\simeq 0.795$, $0.845$, and $1.435$ are of first order, while the rest two at $h\simeq 1.325$ and $2.495$ are of second order. $m_{\mathbf{Q}_{\eta}}$ plotted in the middle panel show that the four phases for $h\lesssim 1.435$ are $4Q$ states, the phase for $1.435\lesssim h\lesssim 2.495$ is a single-$Q$ ($1Q$) conical state, and that for $h\gtrsim 2.495$ is a FFM state. In the $4Q$ states, the amplitudes of $m_{\mathbf{Q}_{\eta}}$ are equal at zero field, while they split into two groups for nonzero fields. We note that the $1Q$ conical phase breaks $C_{2}$ rotational symmetry spontaneously by choosing one of two equivalent wave vectors $\mathbf{Q}_{3}$ and $\mathbf{Q}_{4}$ (we denote the chosen wave vector as $\mathbf{Q}_{4}$ in the figure).

As shown in the bottom panel of Fig. 3(b), $N_{\mathrm{m}}$ is halved through the first-order phase transition at $h\simeq 0.845$ and vanishes through the second-order one at $h\simeq 1.325$. As plotted in the inset, $d_{\rm m}$ gets shorter when approaching $h\simeq 1.325$, similar to the cases of $\mathbf{h}_{001}$ with $h\simeq 0.575$ and $1.395$. This also suggests a topological transition by pair annihilation. On the other hand, $\chi_{\rm sc}$ has a nonzero value in all the $4Q$ states. Notably, $\chi_{\rm sc}$ is almost doubled at $h\simeq 0.845$ where $N_{\rm m}$ is halved, and rapidly decreases through the phase transition at $h\simeq 1.325$ where $N_{\rm m}$ vanishes. We will discuss the relation to the flux flows in Sec. VI.1.

Finally, we discuss the results for the [111] field, $\mathbf{h}_{111}=\frac{1}{\sqrt{3}}(h,h,h)$, shown in Fig. 3(c). $m$, $\chi$, and $m_{\mathbf{Q}_{\eta}}$ in the top and middle panels signal two first-order phase transitions at $h\simeq 0.775$ and $0.845$ among the $4Q$-HLs, a first-order one to the $1Q$ conical state at $h\simeq 1.095$, and a second-order one to a FFM state at $h\simeq 2.595$. In the $4Q$ states, all four $m_{\mathbf{Q}_{\eta}}$ become inequivalent for $0.775\lesssim h\lesssim 0.845$, while two of them have the same amplitudes for $h\lesssim 0.775$ and $0.845\lesssim h\lesssim 1.095$. This indicates that the $4Q$ state for $0.775\lesssim h\lesssim 0.845$ has lower symmetry compared to the other two $4Q$ states, while $C_{3}$ rotational symmetry around the [111] axis ($\parallel\mathbf{Q}_{4}$) is broken in all three phases except at $h=0$.

As plotted in the bottom panel in Fig. 3(c), $N_{\rm m}$ is reduced to half through the first-order phase transition at $h\simeq 0.845$. This leads to the enhancement of $\chi_{\rm sc}$, similar to the case with $\mathbf{h}_{110}$. In the [111] field, however, the system does not exhibit a continuous phase transition that might be ascribed to the topological phase transition. This is presumably due to the fact that the $1Q$ conical state is more stable down to a lower field, compared to the [001] and [110] cases, since the field is applied in parallel to one of the wave vectors, $\mathbf{Q}_{4}$.

V.2 3$Q$ case

Next, we discuss the $3Q$ case with the ordering vectors shown in Fig. 1(b). Figure 4 summarizes the results for the $3Q$-HL at $D=0.3$ and $K=0.7$.

First, we discuss the results for the [001] field, $\mathbf{h}_{001}=(0,0,h)$, shown in Fig. 4(a). As shown in the top and middle panels, $m$, $\chi$, and $m_{\mathbf{Q}_{\eta}}$ signal at least five phase transitions: first-order ones at $h\simeq 0.275$, $0.695$, and $0.775$, and second-order ones at $h\simeq 1.035$ and $2.595$. The four low-field phases for $h\lesssim 1.035$ are $3Q$ states with nonzero three $m_{\mathbf{Q}_{\eta}}$, the phase for $1.035\lesssim h\lesssim 2.595$ is a $1Q$ conical state with only $m_{\mathbf{Q}_{3}}\neq 0$, and that for $h\gtrsim 2.595$ is a FFM state. Furthermore, when we look closer $m_{\mathbf{Q}_{\eta}}$, we find that $m_{\mathbf{Q}_{1}}$ becomes inequivalent to $m_{\mathbf{Q}_{2}}$ at $0.695\lesssim h\lesssim 0.775$, whereas $m_{\mathbf{Q}_{1}}=m_{\mathbf{Q}_{2}}$ in the other three $3Q$ states. These 3$Q$ states are also distinguished by the higher Fourier components of the spin structure factor $S(\mathbf{q})$ and the structure factor of the local scalar spin chirality (see Appendix B).

As shown in the bottom panel of Fig. 4(a), $N_{\rm m}$ is unchanged in the three low-field $3Q$ phases, but it vanishes through the first-order phase transition at $h\simeq 0.775$. On the other hand, $\chi_{\rm sc}$ increases in the two low-field phases, while it rapidly decreases in the third phase and vanishes through the second-order phase transition to the $1Q$ conical state at $h\simeq 1.035$. The change of $\chi_{\rm sc}$ in the $3Q$-HL phases is accounted for by the change in the lengths of the flux flows connecting the monopoles and anti-monopoles, similar to the $4Q$ case in Sec. V.1 (see Sec. VI.2).

Next, we discuss the results for the [110] field, $\mathbf{h}_{110}=\frac{1}{\sqrt{2}}(h,h,0)$, shown in Fig. 4(b). As plotted in the top panel, the data of $m$ and $\chi$ signal seven first-order phase transitions at $h\simeq 0.405$, $0.525$, $0.725$, $0.975$, $0.995$, $1.125$, and $1.245$, and a second-order phase transition at $h\simeq 2.415$. In addition, $m_{\mathbf{Q}_{\eta}}$ in the middle panel and $N_{\rm m}$ in the bottom panel indicate additional phase transitions at $h\simeq 0.445$, $0.745$, $0.945$, and $1.355$. $m_{\mathbf{Q}_{\eta}}$ shows that all the phases for $h\lesssim 1.245$ are $3Q$ states, the two phases for $1.245\lesssim h\lesssim 2.415$ are $2Q$ states, and the phase for $h\simeq 2.415$ is a FFM state. We note that $m_{\mathbf{Q}_{1}}$ becomes inequivalent to $m_{\mathbf{Q}_{2}}$ in the $3Q$ states for $0.405\lesssim h\lesssim 0.525$ and $0.725\lesssim h\lesssim 0.975$, and the $2Q$ state for $1.245\lesssim h\lesssim 1.355$. This indicates spontaneous symmetry breaking by choosing one of the two equivalent wave vectors in these states.

Within the $3Q$ phases for $h\lesssim 1.245$, $N_{\rm m}$ changes in a complicated manner, as plotted in the bottom panel of Fig. 4(b): In contrast to the other cases, $N_{\rm m}$ is not reduced monotonically but changes from 8, 6, 10, 6, to 2 stepwisely. By tracing $d_{\rm m}$ plotted in the insets, we find that the three phase transitions at $h\simeq 0.445$, $0.745$, and $0.945$ appear to be topological ones caused by pair annihilation of monopoles and anti-monopoles ($d_{\rm m}$ does not change from $1$ before the transitions since the lattice spacing is larger than the positional changes of monopoles and anti-monopoles; see Sec. VI.2 for the details). The net scalar spin chirality $\chi_{\rm sc}$ is nonzero in all the $3Q$-HLs for $0<h\lesssim 1.245$. It exhibits a broad peak at $h\sim 1$ and vanishes through the first-order phase transition to the $2Q$ state at $h\simeq 1.245$.

Finally, we discuss the results for the [111] field, $\mathbf{h}_{111}=\frac{1}{\sqrt{3}}(h,h,h)$, shown in Fig. 4(c). In this case, $m$ and $\chi$ plotted in the top panel signal four first-order phase transitions at $h\simeq 0.585$, $1.105$, $1.135$, and $1.255$, and two second-order ones at $h\simeq 1.305$ and $2.335$. In addition, $N_{\rm m}$ in the bottom panel indicates an additional phase transition at $h\simeq 0.495$. $m_{\mathbf{Q}_{\eta}}$ in the middle panel shows that all the phases for $h\lesssim 2.335$ are $3Q$ states, while the phase for $h\gtrsim 2.335$ is a FFM state. All the $3Q$ phases have the equal amplitudes for the three $m_{\mathbf{Q}_{\eta}}$; namely, they retain $C_{3}$ rotational symmetry with respect to the [111] axis.

As plotted in the bottom panel of Fig. 4(c), $N_{\rm m}$ is nonzero in the $3Q$ phases below $h\simeq 1.305$. By monitoring $d_{\rm m}$ plotted in the insets, we find that the transitions at $h\simeq 0.495$ and $1.305$ appears to be topological ones by pair annihilation of monopoles and anti-monopoles. $\chi_{\mathrm{sc}}$ is nonzero for all the $3Q$-HLs but decreases rapidly through the second-order phase transition at $h\simeq 1.305$ where $N_{\mathrm{m}}$ vanishes. We will discuss the details in Sec. VI.2.

VI topological phase transitions by pair annihilation of monopoles and anti-monopoles

In Sec. V, we found several phase transitions in the $4Q$- and $3Q$-HL phases where no discontinuous changes are observed in $m$ and $m_{\mathbf{Q}_{\eta}}$ but $N_{\rm m}$ changes. These suggest continuous phase transitions with a topological change caused by pair annihilation of monopoles and anti-monopoles. Such topological transitions under the [001] field were discussed for an ansatz of the $3Q$-HL state in the continuum limit Zhang et al. (2016) and also for a metastable $3Q$-HL in the model in Eq. (2) Okumura et al. . Our results in Sec. V, however, appear to offer several examples in the ground state for both $4Q$- and $3Q$-HLs. In this section, we analyze these phase transitions by tracing the positions of monopoles and anti-monopoles in real space. In Secs. VI.1 and VI.2, we present the results for the $4Q$ and $3Q$ cases, respectively.

VI.1 4$Q$ case

In Sec. V.1, we found three possible topological phase transitions in the $4Q$-HLs: Two are at $h\simeq 0.575$ and $1.395$ for the [001] field and the other is at $h\simeq 1.325$ for the [110] field. We discuss how the monopoles and anti-monopoles move and pair annihilate as a function of the field strength through each transition.

In the case of the [001] field, $N_{\rm m}$ changes from $16$ to $8$ at $h\simeq 0.575$ and from $8$ to [math] at $h\simeq 1.395$, both suggesting four pairs of monopoles and anti-monopoles annihilate simultaneously at the phase transition. They are visualized in real space in Fig. 5. At zero field, the monopoles and anti-monopoles form two inter-penetrating body-centered-cubic lattices with $d_{\mathrm{m}}=2\sqrt{3}$ as shown in Fig. 5(a). While increasing $h$, half of the monopoles and anti-monopoles move toward each other in the field direction, forming four pairs. When approaching to the critical field, $d_{\rm m}$ given by the four pairs is reduced to $1$ as shown in Fig. 5(b) at $h=0.57$, and then, becomes [math], which is the pair annihilation at the critical field $h\simeq 0.575$. In the higher-field region, the remaining monopoles and anti-monopoles are paired along the field direction again, as exemplified in Fig. 5(c) at $h=0.60$. In this case, $d_{\rm m}$ is reduced from $3$ to $1$ as shown in Fig. 5(d) at $h=1.39$, and finally, becomes [math] at the critical field $h\simeq 1.395$ by the pair annihilation. The changes of $d_{\rm m}$ were plotted in the inset of the bottom panel in Fig. 3(a).

The movement of the monopoles and anti-monopoles explains the behavior of $\chi_{\mathrm{sc}}$ plotted in the bottom panel of Fig. 3(a). When approaching the topological transition at $h\simeq 0.575$ by increasing $h$, $\chi_{\mathrm{sc}}$ decreases [$-\chi_{\mathrm{sc}}$ increases in Fig. 3(a)]. This is understood by the decrease of $d_{\mathrm{m}}$ with the flux flows in the same direction of the magnetic field: The decrease of $d_{\mathrm{m}}$ reduces the positive contribution to $\chi_{\mathrm{sc}}$, which leads to the net decrease in $\chi_{\mathrm{sc}}$. On the other hand, $\chi_{\mathrm{sc}}$ increases ($-\chi_{\mathrm{sc}}$ decreases) near the other topological transition at $h\simeq 1.395$. This is due to the decrease of $d_{\mathrm{m}}$ with the flux flows in the opposite direction to the magnetic field.

Similarly, $N_{\rm m}$ changes from $8$ to [math] for the [110] field through the phase transition at $h\simeq 1.325$. The change of the positions of monopoles and anti-monopoles is shown in Fig. 6, where $d_{\rm m}$ changes in a similar manner to the case of the [001] field at $h\simeq 1.39$ in Figs. 5(c) and 5(d); see also the inset of the bottom panel of Fig. 3(b). The only difference from the [001] case is in the direction of collisions. The corresponding reduction of the lengths of the flux flows is also related to the suppression of $\chi_{\mathrm{sc}}$ in Fig. 3(b) since the fluxes $\bm{\Omega}_{p}$ have the positive component in the opposite direction to the field.

VI.2 3$Q$ case

In the case of the $3Q$-HLs, we identified totally five possible topological transitions in Sec. V.2. Three of them are at $h\simeq 0.445$, $0.745$, and $0.945$ for the [110] field, and the rest two are at $h\simeq 0.495$ and $1.305$ for the [111] field. Figure 7 shows the real-space pictures for the [110] field. In a low field, there are totally eight monopoles and anti-monopoles as shown in Fig. 7(a) for $h=0.41$, but one pair annihilates through the transition at $h\simeq 0.445$ as shown in Figs. 7(b) and 7(c). We note that the directions of pairs change within the same $3Q$ state with spontaneous symmetry breaking with respect to the [100] and [010] directions. In the next topological transition at $h\simeq 0.745$, $N_{\mathrm{m}}$ changes from $10$ to $6$, where two pairs of monopoles and anti-monopoles annihilate as shown in Figs. 7(d) and 7(e). Through these transitions, $d_{\rm m}$ does not change from $1$, since the distance for the pairs that survive is already $1$ before the transition. Finally, other two pairs annihilate and $N_{\mathrm{m}}$ is reduced to 2 at $h\simeq 0.945$ as shown in Fig. 7(f). See also the insets in the bottom panel of Fig. 4(b).

Finally, we present the results for the [111] field in Fig. 8. In this case, there are four pairs of monopoles and anti-monopoles in the low-field phase, and their distance gets shorter as demonstrated in Figs. 8(a) and 8(b). Through the topological transition at $h\simeq 0.495$, three of four annihilate as shown in Fig. 8(c). Finally, the remaining pair gets closer and pair annihilate through the transition at $h\simeq 1.305$ as shown in Figs. 8(d) and 8(e). The rapid decrease of $d_{\mathrm{m}}$ explains the rapid suppression of $\chi_{\mathrm{sc}}$ while approaching the topological phase transition at $h\simeq 1.305$ in Fig. 4(c).

VII Concluding remarks

In conclusion, we have investigated the magnetic HLs in the effective spin model with long-range interactions reflecting the itinerant nature of electrons. We found that both $4Q$- and $3Q$-HLs are stabilized even at zero magnetic field by the synergy between the DM-type interactions from the spin-orbit coupling and the multiple-spin interactions from the spin-charge coupling. The results are in stark contrast to the previous studies for the localized spin models with short-range interactions, where the HLs are stable only in a field. Furthermore, our HLs may have much shorter periods compared to the previous ones; the periods in our HLs are dictated by nesting properties of the Fermi surface, whereas those in the previous studies are given by the competition between the ferromagnetic exchange interaction and the DM interaction. We also clarified the effect of an external magnetic field on the HLs. We showed that both $4Q$ and $3Q$ cases exhibit a variety of successive phase transitions depending on the field direction, including the transitions to $2Q$ and $1Q$ states. Interestingly, among them, we found several topological phase transitions where the number of monopoles and anti-monopoles changes by the pair annihilation. We explicitly showed how the pair annihilation occurs by tracing the real-space positions of monopoles and anti-monopoles on the discrete lattice.

As mentioned in the introduction, $3Q$- and $4Q$-HLs were recently discovered in MnSi1-xGex Tanigaki et al. (2015); Kanazawa et al. (2017); Fujishiro et al. (2019). They are stable even in the absence of the magnetic field and have much shorter periods compared to the conventional skyrmion lattices, for instance, in MnSi, and evaded the understanding from the conventional spin models with short-range two-spin interactions. A scenario was recently proposed based on short-range four-spin and six-spin interactions including the scalar spin chirality Grytsiuk et al. (2019). Our finding suggests another scenario by emphasizing the important role of itinerant nature of electrons. To test our scenario, it is necessary to clarify the electronic structure in the real compounds, e.g., by the angle-resolved photoemission spectroscopy and the de Haas-van Alphen effect. First-principles calculations would also be helpful, while it is not straightforward to precisely predict the relevant wave numbers in the complicated multiorbital systems with electron correlations, in particular, chemically doped materials like MnSi1-xGex. It would also be interesting to test our scenario for the short-period skyrmion lattice recently discovered in EuPtSi Kakihana et al. (2018); Kaneko et al. (2019); Takeuchi et al. (2019). We note that a similar scenario (without the DM-type interaction) was recently discussed for the swirling spin textures in a centrosymmetric triangular magnet Gd2PdSi3 Kurumaji et al. (2019).

On the other hand, in the magnetic field, our results suggest that the $4Q$ and $3Q$ states exhibit a nonzero topological Hall effect through the nonzero scalar spin chirality $\chi_{\mathrm{sc}}$. Our results also indicate that $\chi_{\rm sc}$ changes drastically corresponding to the modulation of the magnetic textures including the topological transitions by pair annihilations of monopoles and anti-monopoles. Experimentally, interesting behaviors were observed in a wide range of field and temperature, even with the sign change of the topological Hall resistivity Fujishiro et al. (2019). Assuming our scenario based on the itinerant nature of electrons, it will be important to take into account the realistic electronic band structures in the magnetic field for detailed comparison between theory and experiment. In particular, it is worth studying how the modulations of the Fermi surfaces and corresponding $\mathbf{Q}_{\eta}$ modify the phase diagrams in the magnetic field. Moreover, thermal fluctuations might also play an important role. We leave the finite-temperature study as a future work, as it requires sophisticated Monte Carlo simulations beyond the simulated annealing to resolve competing phases.

Acknowledgements.

We would like to thank N. Kanazawa and K. Shimizu for fruitful discussions. This research was supported by Grants-in-Aid for Scientific Research under Grants No. JP19H05825 and No. JP18K13488, JST CREST (JP- MJCR18T2), and the Chirality Research Center in Hiroshima University and JSPS Core-to-Core Program, Advanced Research Networks. S. O. was supported by JSPS through the research fellowship for young scientists.

Appendix A Difference among the 4$Q$ states

In this section, we discuss the difference among the $4Q$ states found in Sec. V.1. In the case of the [001] field in Fig. 3(a), we found four $4Q$ states below $h\simeq 2.335$, all of which have the equal amplitudes of the four $m_{\mathbf{Q}_{\eta}}$. Two of them are distinguished by the number of monopoles and anti-monopoles, $N_{\rm m}$: the $4Q$-HL state with $N_{\rm m}=16$ below $h\simeq 0.575$ and the $4Q$ state with $N_{\rm m}=0$ above $h\simeq 1.395$. The rest two have the same $N_{\rm m}=8$, but we find that they show different values in higher harmonics in the spin structure factor $S(\mathbf{q})$ in Eq. (9): By calculating $S(\mathbb{Q}_{\eta}+\mathbb{Q}_{\eta^{\prime}})$ where $\eta=1,2,3,4$ and $\eta\neq\eta^{\prime}$, we find that $S(\mathbb{Q}_{2}+\mathbb{Q}_{3})$ and $S(\mathbb{Q}_{3}+\mathbb{Q}_{1})$ have the equal amplitudes for $0.575\lesssim h\lesssim 0.595$, but $S(\mathbb{Q}_{2}+\mathbb{Q}_{3})\neq S(\mathbb{Q}_{3}+\mathbb{Q}_{1})$ for $0.595\lesssim h\lesssim 1.395$. The differences among the four $4Q$ states are summarized in Fig. 9(a).

In the case of the [110] field in Fig. 3(b), we also find four $4Q$ states below $h\simeq 1.435$ that share the same symmetry in terms of $m_{\mathbf{Q}_{\eta}}$. In this case again, $N_{\rm m}$ distinguishes two of them: the $4Q$-HL state with $N_{\rm m}=8$ for $0.845\lesssim h\lesssim 1.325$ and the $4Q$ state with $N_{\rm m}=0$ for $1.325\lesssim h\lesssim 1.435$. In order to distinguish the rest two for $h\lesssim 0.845$, we calculate the spin structure factor with the spin component parallel to the [110] field defined by

[TABLE]

where $S^{\mu\nu}(\mathbf{q})$ is the matrix form of the spin structure factor defined by

[TABLE]

We find that the higher harmonics along the field direction, $S^{\parallel}(2\mathbf{Q}_{1}+2\mathbf{Q}_{2})$, is nonzero in the $4Q$ state for $0.795\lesssim h\lesssim 0.845$, while it vanishes for $h\lesssim 0.795$. The differences among the four $4Q$ states are summarized in Fig. 9(b).

Finally, in the case of the [111] field in Fig. 3(c), we find three $4Q$ states. In this case, $m_{\mathbf{Q}_{\eta}}$ distinguishes the intermediate phase for $0.845\lesssim h\lesssim 1.095$, as mentioned in Sec. V.1. The rest two can be distinguished by $N_{\rm m}$. See Fig. 9(c).

Appendix B Difference among the 3$Q$ states

In this section, we discuss the difference among the $3Q$ states found in Sec. V.2. In the case of the [001] field in Fig. 4(a), we find four $3Q$ states. One of them for $0.775\lesssim h\lesssim 1.035$ is distinguished from the others by $N_{\rm m}$: the $3Q$ state with $N_{\rm m}=0$ above $h\simeq 0.775$. In addition, as mentioned in Sec. V.2, $m_{\mathbf{Q}_{\eta}}$ distinguishes the $3Q$-HL state for $0.695\lesssim h\lesssim 0.775$. We find a difference between the rest two in higher harmonics in the structure factor of the scalar spin chirality in Eq. (15). Specifically, we calculate the component perpendicular to the [001] field defined by

[TABLE]

where $\chi_{\mathrm{sc}}^{\mu\nu}(\mathbf{q})$ is the matrix form of the structure factor defined by

[TABLE]

We find that $\chi_{\mathrm{sc}}^{\perp}(2\mathbf{Q}_{1})$ and $\chi_{\mathrm{sc}}^{\perp}(2\mathbf{Q}_{2})$ are nonzero and have the equal amplitudes in the $3Q$-HL state for $0.275\lesssim h\lesssim 0.695$, but vanish below $h\simeq 0.275$. The differences among the four $3Q$ states are summarized in Fig. 10(a).

Next, in the case of the [110] field in Fig. 4(b), we found ten $3Q$ states below $h\simeq 1.245$ with a variety of $N_{\rm m}$. The $3Q$-HL state with $N_{\rm m}=10$ for $0.725\lesssim h\lesssim 0.745$ is distinguished from others, but $N_{\rm m}=8$ for the two states below $h\simeq 0.445$, $N_{\rm m}=6$ for the three states for $0.445\lesssim h\lesssim 0.725$ and $0.745\lesssim h\lesssim 0.945$, and $N_{\rm m}=2$ for the rest four for $0.945\lesssim h\lesssim 1.245$. As mentioned in Sec. V.2, the two with $N_{\rm m}=8$ below $h\simeq 0.445$ and the two with $N_{\rm m}=6$ for $0.445\lesssim h\lesssim 0.725$ are distinguished by $m_{\mathbf{Q}_{\eta}}$. Similarly, the two with $N_{\rm m}=2$ for $0.945\lesssim h\lesssim 0.975$ and $1.125\lesssim h\lesssim 1.245$ are distinguished from other two by $m_{\mathbf{Q}_{\eta}}$. See Fig. 10(b).

In order to distinguish the rest, we calculate higher harmonics in the spin structure factor similar to the $4Q$ case in Appendix A. For the $3Q$-HL states with $N_{\rm m}=6$, $S^{\parallel}(\mathbb{Q}_{2}+\mathbb{Q}_{3})$ and $S^{\parallel}(\mathbb{Q}_{3}+\mathbb{Q}_{1})$ have the equal amplitudes for $0.745\lesssim h\lesssim 0.945$, but do not for $0.445\lesssim h\lesssim 0.525$. Meanwhile, for the two states with $N_{\rm m}=2$ and $m_{\mathbf{Q}_{1}}\neq m_{\mathbf{Q}_{2}}$, $S^{\parallel}(2\mathbb{Q}_{2}+2\mathbb{Q}_{3})=S^{\parallel}(2\mathbb{Q}_{3}+2\mathbb{Q}_{1})$ for $1.125\lesssim h\lesssim 1.245$ while $S^{\parallel}(2\mathbb{Q}_{2}+2\mathbb{Q}_{3})\neq S^{\parallel}(2\mathbb{Q}_{3}+2\mathbb{Q}_{1})$ for $0.945\lesssim h\lesssim 0.975$. Furthermore, for the two states with $N_{\rm m}=2$ and $m_{\mathbf{Q}_{1}}=m_{\mathbf{Q}_{2}}$ the component perpendicular to the [110] field, $S^{zz}(2\mathbb{Q}_{1}+2\mathbb{Q}_{2})$, is zero for $0.975\lesssim h\lesssim 0.995$, but nonzero for $0.995\lesssim h\lesssim 1.125$. All the differences among the ten $3Q$ states are summarized in Fig. 10(b).

Finally, in the case of the [111] field in Fig. 4(c), we found seven $3Q$ states below $h\simeq 2.335$, all of which have the equal amplitudes of the three $m_{\mathbf{Q}_{\eta}}$. Two of them are distinguished by $N_{\rm m}$: the $3Q$-HL state with $N_{\rm m}=8$ below $h\simeq 0.495$ and the $3Q$ state with $N_{\rm m}=0$ above $h\simeq 1.305$. The rest five have the same $N_{\rm m}=2$, but two of them for $1.105\lesssim h\lesssim 1.135$ and $1.255\lesssim h\lesssim 1.305$ show nonzero values in higher harmonics in the spin structure factor $S^{\perp}(\mathbf{Q}_{1}+\mathbf{Q}_{2}+\mathbf{Q}_{3})$, which is the component perpendicular to the [111] field given by

[TABLE]

Furthermore, $S^{\perp}(\mathbf{Q}_{1}+\mathbf{Q}_{2})=S^{\perp}(\mathbf{Q}_{2}+\mathbf{Q}_{3})$ for $1.105\lesssim h\lesssim 1.135$, but $S^{\perp}(\mathbf{Q}_{1}+\mathbf{Q}_{2})\neq S^{\perp}(\mathbf{Q}_{2}+\mathbf{Q}_{3})$ for $1.255\lesssim h\lesssim 1.305$. The rest three $3Q$-HLs (for $0.495\lesssim h\lesssim 0.585$, $0.585\lesssim h\lesssim 1.105$, and $1.135\lesssim h\lesssim 1.255$) with $N_{\rm m}=2$ cannot be distinguished within the present analyses although further higher harmonics may tell the difference. See Fig. 10(c).