Theory of Collective Excitations in the Quadruple-$Q$ Magnetic Hedgehog Lattices

Abstract

Hedgehog and antihedgehog spin textures in magnets behave as emergent monopoles and antimonopoles, which give rise to astonishing transport and electromagnetic phenomena. Using the Kondo-lattice model in three dimensions, we theoretically study collective spin-wave excitation modes of magnetic hedgehog lattices which have recently been discovered in itinerant magnets such as MnSi_1-xGe_x and SrFeO₃. It is revealed that the spin-wave modes, which appear in the sub-terahertz regime, have dominant amplitudes localized at Dirac strings connecting hedgehog-antihedgehog pairs and are characterized by their translational oscillations. It is found that their spectral features sensitively depend on the number and configuration of the Dirac strings and, thus, can be exploited for identifying the topological phase transitions associated with the monopole-antimonopole pair annihilations.

Figure 1: (a), (b) Schematics of hedgehog (a) and antihedgehog (b) spin textures. Green arrows represent emergent magnetic fields around the point defects (Bloch points), and arrows on the cube corners represent the localized spins. (c) Spin configuration of String A. (d), (e) Spatial distributions of the Bloch points of hedgehogs (magenta) and antihedgehogs (cyan) in the 4$Q$-HLs with (d) $N_{\rm m}$=16 and (e) $N_{\rm m}$=8. Red and blue lines represent the Dirac strings with vorticities of $-1$ and $+1$. (f) Four propagation vectors of magnetic helices constituting the 4$Q$-HLs. (g) Magnetic helix with a propagation vector $\bm{Q}_{\nu}$ where $\theta_{\nu}$ represents the phase degree of freedom with respect to the translation of helix. (h) Fermi surface for the kinetic term of $\mathcal{H}_{\rm KLM}$ in Eq. (1). (i) Two-dimensional cross section of the Fermi surface with nesting vectors $\bm{Q}_{1}$ and $\bm{Q}_{2}$.

Magnetic monopoles, elementary particles with isolated magnetic charges in three dimensions, have attracted continuous research interest since Dirac’s first proposal in 1931 [1]. Although real magnetic monopoles have not been observed in nature, scientists have discovered quasiparticles that virtually behave as magnetic monopoles in condensed-matter systems such as Dirac-electron materials [2] and spin-ice pyrochlores [3, 4]. Interesting physical phenomena associated with such emergent magnetic monopoles, e.g., fractional excitations [5, 3, 6, 7], anomalous Hall effects [8, 9], and anomalous Nernst effects [10, 11], have been discussed and/or observed in these systems.

Recently, another condensed-matter system that realizes emergent magnetic monopoles was discovered [12, 13, 14]. It was revealed that an itinerant chiral magnet MnGe hosts a periodic array of magnetic hedgehogs and antihedgehogs called magnetic hedgehog lattice [15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29]. These hedgehogs and antihedgehogs can be regarded as magnetic monopoles and antimonopoles, respectively, as they behave as sources and sinks of emergent magnetic fields acting on the conduction electrons via exchange coupling through the Berry-phase mechanism [Figs. 1(a)-(c)].

The magnetic hedgehog lattice in MnGe is described by a superposition of spin helices with cubic three propagation vectors and, thereby, is referred to as a triple-$Q$ hedgehog lattice (3$Q$-HL). Recently, it was experimentally discovered that substitution of Ge with Si transforms this 3$Q$-HL into another type of hedgehog lattice [Figs. 1(d),(e)]. A small-angle neutron-scattering experiment revealed that the hedgehog lattice in MnSi_1-xGe_x is characterized by a superposition of spin helices with tetrahedral four propagation vectors [Fig. 1(f)], i.e., the quadruple-$Q$ hedgehog lattice (4$Q$-HL) [30]. The 4$Q$-HL was observed also in the perovskite ferrite SrFeO₃ [31].

In fact, there are several kinds of 4$Q$-HL states with different number of hedgehogs and antihedgehogs in the magnetic unit cell because the spin structure of 4$Q$-HLs are characterized not only by the propagation vectors $\bm{Q}_{\nu}$ ($\nu$=1-4) [32, 33, 34] but also by other internal degrees of freedom, e.g., net magnetization [32], chirality [35], and relative phases $\theta_{\nu}$ of the superposed spin helices [Fig. 1(g)] [36, 37]. Therefore, we expect that MnSi_1-xGe_x and SrFeO₃ offer precious opportunities to study the tunability of transport and thermoelectric properties via the field-induced variation of topological nature of the hedgehog lattices [36, 38, 39].

The dynamical phenomena associated with these emergent magnetic monopoles are more interesting. From the electromagnetic duality in Maxwell’s equations, the dynamics of magnetic monopoles is expected to induce an electric field, as the dynamics of electric charge induces a magnetic field. Through inducing the electric field while moving, the emergent monopoles behaves as dyons, i.e., hypothetical particles in the grand unified theory [40]. Novel transport phenomena and optical/microwave device functions are expected to emerge from such emergent dyons [41, 42, 43, 44, 45]. In order to study new physical phenomena and device functions originating from the dynamics of such emergent monopoles and dyons, knowledge of their intrinsic excitations is essentially important.

In this Letter, we theoretically study collective excitation modes of the 4$Q$-HLs in itinerant chiral magnets using a microscopic spin-charge coupled model in three dimensions. We construct a model on the chiral-cubic lattice so as to reproduce the distinct two types of 4$Q$-HLs with very short periods of a few atomic sites, which are realized in MnSi_1-xGe_x and SrFeO₃, based on the microscopic insights into these materials. We first reveal that the nesting of Fermi surfaces can work as a principal mechanism for stabilizing the 4$Q$-HL textures. Then we find that the 4$Q$-HLs have characteristic collective excitation modes associated with translational motion of Dirac strings. Because these collective excitation modes are dominated by the Dirac strings, presence or absence of the modes sensitively depends on the number of monopole-antimonopole pairs and their spatial configuration. Indeed, we demonstrate that a certain mode disappears upon a topological phase transition associated with vanishing of Dirac strings due to the field-induced monopole-antimonopole pair annihilations. This finding will contribute to characterization of the magnetic topologies in hedgehog-hosting magnets and will be a basis for studying the monopole-induced physical phenomena in matters.

We start with the Kondo-lattice model on a cubic lattice with the Dzyaloshinskii-Moriya (DM) interactions, which describes the coupling between the conduction electrons and the localized spins in itinerant chiral magnets, and some additional terms. The Hamiltonian is given by $\mathcal{H}=\mathcal{H}_{\rm KLM}+\mathcal{H}_{\rm Zeeman}+\mathcal{H}_{\rm local}$ with

$$\mathcal{H}_{\rm KLM}=-\sum_{i,j,\sigma}t_{ij}\hat{c}_{i\sigma}^{\dagger}\hat{c}_{j\sigma}-J_{\rm K}\sum_{i}\hat{\bm{s}}_{i}\cdot{\bm{S}}_{i},$$ (1)

$$\mathcal{H}_{\rm Zeeman}=-\left[\bm{H}+\bm{h}(t)\right]\cdot\sum_{i}\bm{S}_{i},$$ (2)

$$\mathcal{H}_{\rm local}=\sum_{<i,j>}J_{\rm AFM}\bm{S}_{i}\cdot\bm{S}_{j}-D\sum_{<i,j>}\bm{e}_{ij}\cdot(\bm{S}_{i}\times\bm{S}_{j}),$$ (3)

where $\hat{c}_{i\sigma}^{\dagger}$ ($\hat{c}_{i\sigma}$) denotes a creation (annihilation) operator of a conduction electron with spin $\sigma(=\uparrow,\downarrow)$ at site $i$. The first term of $\mathcal{H}_{\rm KLM}$ describes the kinetic energy of the conduction electrons with the nearest-neighbor hopping $t_{1}(=1)$ and the fourth-neighbor hopping $t_{4}(=-1)$. The second term of $\mathcal{H}_{\rm KLM}$ describes the Kondo exchange coupling between the localized classical spins $\bm{S}_{i}$ $(|\bm{S}_{i}|=1)$ and the conduction-electron spins $\hat{\bm{s}}_{i}=\left(1/2\right)\sum_{\sigma\sigma^{\prime}}\hat{c}_{i\sigma}^{\dagger}{\bm{\sigma}}_{\sigma\sigma^{\prime}}\hat{c}_{i\sigma^{\prime}}$ where $\bm{\sigma}$ represents the vector of Pauli matrices. The term $\mathcal{H}_{\rm Zeeman}$ denotes Zeeman coupling associated with the external magnetic field, while the term $\mathcal{H}_{\rm Local}$ describes the antiferromagnetic (AFM) exchange interactions and the DM interactions between the nearest-neighbor localized spins, where ${\bm{e}}_{ij}$ represents the bond-directional unit vectors. These interactions originate from the hybridizations among the localized orbitals and the spin-orbit coupling [46]. We set $J_{\rm K}$=0.8 and $J_{\rm AFM}$=0.0008, and we examine both the cases with and without DM interactions, i.e., $D=0$ and $D=0.0002$.

The 4$Q$-HL texture is described by a superposition of four magnetic helices governed by the nested Fermi surfaces of $\mathcal{H}_{\rm KLM}$ [Figs. 1(h) and (i)]. More specifically, the propagation vectors of helices $\{\bm{Q}_{\nu}\}$ coincide the nesting vectors [Fig. 1(i)]. This is because the effective two-body exchange interactions between the localized spins are governed by the bare susceptibility $\chi^{0}_{\bm{q}}$, which has maxima at momenta $\bm{q}$ that correspond to the nesting wavenumbers [47]. Indeed, when $J_{\rm K}$ is weak, the second-order perturbation theory gives $-(J_{\rm K}/2)^{2}\sum_{{\bm{q}}\in{\rm BZ}}\chi^{0}_{\bm{q}}|\bm{S}_{\bm{q}}|^{2}$ as the effective interaction.

The tetrahedral four propagations vectors $\{\bm{Q}_{\nu}\}$ are represented by $\bm{Q}_{1}=(Q_{\rm abs},Q_{\rm abs},Q_{\rm abs})$, $\bm{Q}_{2}=(-Q_{\rm abs},-Q_{\rm abs},Q_{\rm abs})$, $\bm{Q}_{3}=(-Q_{\rm abs},Q_{\rm abs},-Q_{\rm abs})$, and $\bm{Q}_{4}=(Q_{\rm abs},-Q_{\rm abs},-Q_{\rm abs})$. In the present study, we set the chemical potential $\mu=-3.79$, which gives $Q_{\rm abs}\approx\pi/4$. This wavenumber corresponds to a spatial period of $\lambda=2\pi a/\sqrt{3}Q_{\rm abs}\sim$2.15 nm, if we assume the lattice constant of $a=0.465$ nm, which reproduces the experimentally observed $\lambda=$1.9-2.1 nm for MnSi_1-xGe_x [30, 27]. Note that the following results are robust, and a fine tuning of the chemical potential is not required [48].

Figure 2: (a) Ground-state phase diagram as a function of the external magnetic field $H_{z}$, which contains the 4$Q$-HL with $N_{\rm m}$=16, another 4$Q$-HL with $N_{\rm m}$=8, the intermediate 4$Q$ (I-4$Q$) phase, and the forced ferromagnetic phase. (b) Imaginary part of the longitudinal dynamical magnetic susceptibility for various values of $H_{z}$. The DM parameter is set to be $D=0.0002$.

We investigate the ground-state phase diagram of the Hamiltonian $\mathcal{H}$ by using the adiabatic Landau-Lifshitz-Gilbert equation, $\dot{\bm{S}}_{i}=-{\bm{S}}_{i}\times{\bm{H}}_{i}^{\rm eff}+\alpha_{\rm G}{\bm{S}}_{i}\times\dot{\bm{S}}_{i}$, where $\alpha_{\rm G}$ is the Gilbert-damping coefficient. The effective local magnetic field is calculated by the first-order differential of thermodynamic potential as $\bm{H}_{i}^{\rm eff}=-\partial\Omega/\partial\bm{S}_{i}$. For this formulation, we assume an adiabatic condition that dynamics of the conduction electrons are sufficiently fast compared to that of the localized spins and follow them smoothly, so that the conduction electrons are always in equilibrium. The calculations are numerically performed by using the kernel polynomial method combined with the automatic differentiation technique based on the chain rule of Chebyshev polynomials [51, 52, 53, 54, 55]. This method is known to be powerful when simulating the low-energy dynamics in spin-charge coupled systems [56, 55, 57].

The phase diagram obtained for a system of $16^{3}$ sites is shown in Fig. 2(a). We also confirm that the phase diagram is not changed in an analysis of a larger size system with $32^{3}$ sites. Here the DM parameter is set to be $D=0.0002$. When $H_{z}$=0, a 4$Q$-HL state with $N_{\rm m}=16$ appears where $N_{\rm m}$ is the total number of hedgehogs and antihedgehogs in the magnetic unit cell. It is worth mentioning that the stable 4$Q$-HL state depends on the presence or absence of the DM interaction [32, 35]. When the DM interaction is absent ($D$=0), the 4$Q$-HL consisting of two right-handed helices and two left-handed helices is stabilized. On the contrary, the 4$Q$-HL consisting of four-righthanded helices is stabilized in the presence of the DM interaction ($D=0.0002$). In the following, we focus on the 4$Q$-HL in the latter case. This chiral 4$Q$-HL state contains two inequivalent Dirac strings (A and B), each of which connects a hedgehog and an antihedgehog along the $z$ axis. There are eight strings in total, four each for A and B. Note that the magnitude of DM interaction required for stabilizing the chiral 4$Q$-HL is as small as $D=0.0002$, which is approximately 35% of $J_{1}^{\rm eff}$ and 25% of $J_{\rm AFM}$. Here $J_{1}^{\rm eff}=-(J_{\rm K}/2)^{2}(1/N_{\bm{q}})\sum_{{\bm{q}}\in{\rm BZ}}\chi^{0}_{\bm{q}}e^{i{\bm{q}}\cdot{\bm{e}}_{\gamma}}$ $(\gamma=x,y,z)$ is the coupling constant of the effective nearest-neighbor ferromagnetic exchange interactions mediated by the conduction electrons.

When the external magnetic field is absent ($H_{z}=0$), the Dirac strings A and B have the same length $d_{\rm A}=d_{\rm B}=4a$. As $H_{z}$ is increased, both $d_{\rm A}$ and $d_{\rm B}$ decrease, keeping the relationship $d_{\rm A}\geq d_{\rm B}$. The length is defined by $d=|\bm{R}_{\rm H}-\bm{R}_{\rm AH}|$. Here $\bm{R}_{\rm H}$ ($\bm{R}_{\rm AH}$) is a center position of the unit cube on which the (anti)hedgehog is defined, which takes integer values with the lattice constant $a$ as the unit. When the field strength reaches $H_{z}\approx$0.0046, the length $d_{\rm B}$ reaches zero first, which results in the pair annihilations of the hedgehog and antihedgehog connected by the String B and vanishing of the String B. Consequently, a topological phase transition occurs from the 4$Q$-HL state with $N_{\rm m}$=16 to that with $N_{\rm m}$=8. This phase transition is characterized by a change of the phase $\Theta$, which is defined by $\Theta=\pi-\left|\pi-{\rm mod}\left[\Theta^{\prime},2\pi\right]\right|$ with $\Theta^{\prime}=\sum_{\nu=1}^{4}\theta_{\nu}$, where $\theta_{\nu}$ is the relative phase shift of helix with the propagation vector $\bm{Q}_{\nu}$ [Fig. 1(g)]. With increasing $H_{z}$, $\Theta$ decreases from $\Theta=\pi/3$ and is strongly suppressed around the phase boundary. While $\Theta$ is still finite at the phase boundary, it becomes zero when the system goes a little inside the $N_{\rm m}$=8 phase [Fig. 1(e)] [36]. The ellipticity $\varepsilon_{\nu}$ also exhibits a characteristic change with a minimum around the phase boundary [48].

Table 1: Properties of the collective excitation modes in the 4$Q$-HL states. The second column shows the Dirac strings at which the mode has large amplitude, while the third column shows the relevant 4$Q$-HL state.

Mode	Relevant Dirac strings	Relevant 4$Q$-HL state(s)
L1	-	$N_{\rm m}=16$, $N_{\rm m}=8$
L2	Strings B	$N_{\rm m}=16$
L3	Strings A	$N_{\rm m}=16$, $N_{\rm m}=8$

Now we study the collective excitation modes in these chiral 4$Q$-HL phases by calculating the longitudinal dynamical magnetic susceptibility $\chi_{z}(\omega)$. Time profiles of the total magnetization $S_{z}(t)$ are calculated using the aLLG equation after applying a short-time pulse of $\bm{h}(t)=h_{z}\delta(t)\bm{e}_{z}$ and perform the Fourier transformation of $\Delta S_{z}(t)=S_{z}(t)-S_{z}(0)$ to obtain $\chi_{z}(\omega)$. The Gilbert-damping coefficient is fixed at $\alpha_{\rm G}$=0.04 for the simulations. The obtained spectra of Im$\chi_{z}(\omega)$ show that three intrinsic excitation modes (L1, L2, and L3) appear in the $N_{\rm m}=16$ phase, whereas the mode L2 disappears when the system enters the $N_{\rm m}=8$ phase. Note that these modes appear in the sub-terahertz regime because $\omega$=0.004 in Fig. 2(b) corresponds to 1 THz approximately when we assume $t_{1}=1$ eV.

Figure 3: (a)-(c) Spatial maps of the oscillation amplitudes projected onto the $xy$, $yz$ and $zx$ planes for the L1 mode in the 4$Q$-HL phase with $N_{\rm m}=16$ when $H_{z}=0.001$ and $D=0.0002$. (d)-(f) Those for the L2 mode. (g)-(i) Those for the L3 modes. The red (blue) solid lines represent Strings A (B), while the red (blue) dashed lines connect Strings A (B) to show their spatial configuration seen along the $c$ axis, which are also shown in Figs. 1(d) and 1(e).

A notable property of these longitudinal modes is their spatial localization. Figures 3(a)-(i) show the spatial distribution of oscillation amplitudes for the L1, L2, and L3 modes in the $N_{\rm m}=16$ phase at $H_{z}=0.001$ excited by a time-dependent field $\bm{h}(t)=h_{z}\sin(\omega_{\rm ac}t)\bm{e}_{z}$ with a small amplitude of $h_{z}=10^{-5}$. Here the frequency $\omega_{\rm ac}$ is fixed at the eigenfrequency of the corresponding mode, and we used $\alpha_{\rm G}=0.008$ for the simulations. Similar plots as Fig. 3 for the modes in the $N_{\rm m}=16$ phase at $H_{z}=0$ and those for the modes in the $N_{\rm m}=8$ phase at $H_{z}=0.005$ are presented in Figs. S4 and S5 in the Supplemental Materials [48], respectively. Here sizes of the green balls represent norms of the oscillation amplitudes $\left|\delta{\bm{S}}_{i}\right|=\left[\sum_{\mu=x,y,z}\left(\delta S_{i\mu}\right)^{2}\right]^{1/2}$ during a period of the oscillation. The ball sizes are normalized with respect to the largest value in the unit cell after subtracting a spatially uniform component for clear visibility.

In these figures, we find that the L2 mode is localized at the Strings B, whereas the L3 mode is at the Strings A. Consequently, the L2 mode disappears when the system enters the $N_{\rm m}=8$ phase because the Strings B disappear as the hedgehog-antihedgehog pair annihilations occur upon this phase transition [Fig. 2(b)]. The disappearance of L2 mode manifests the change in magnetic topology from $N_{\rm m}=16$ to $N_{\rm m}=8$. This conversely indicates that this topological transition can be observed in the magnetic resonance spectra. On the contrary, the L3 mode survives after this phase transition [Fig. 2(b)]. With further increasing $H_{z}$, the system enters the intermediate 4$Q$ (I-4$Q$) phase and subsequently the forced ferromagnetic (FFM) phase. In the I-4$Q$ phase, the hedgehog and antihedgehog connected by a Dirac string collide to be merged with their cores sharing the same cubic unit cell [48]. In this phase, there still exists nonzero scalar spin chirality as a residue of Strings A, although it is no longer quantized, and the L3 mode survive because of the remnant of topology. The mode disappears upon the transition to the collinear FFM phase at $H_{z}\approx 0.0081$, in which the scalar spin chirality vanishes [48].

In contrast to the L2 and L3 modes, the L1 mode is not localized, but its oscillation amplitude is widely distributed in the magnetic unit cell. We also mention two aspects of the spin-wave modes. First, the L1 mode is not $C_{4z}$-symmetric with respect to the spatial distribution of the oscillation amplitudes as seen in Figs. 3(b) and 3(c) because of the reduced symmetry of the spin texture [48]. Second, the three spin-wave modes are gapped when $H_{z}$=0 and distinct from those in the 3$Q$ hedgehog lattices discussed in Ref. [41], which are gapless at zero field.

Figure 4: (a),(b) Relationship of the positions of Bloch points in the crystallographic unit cell and local solid angles composed of localized spins. The length of the green arrows represent the magnitude of solid angle spanned by four localized spins, which is roughly proportional to the local emergent magnetic fields. (c) Time profiles of the displacements $\delta_{\rm(A)H}$ of the Bloch points for the String B in the L2 mode. (d) Time profiles of $\delta_{\rm(A)H}$ for the String A in the L3 mode. (e),(f) Schematics of the oscillatory translational motion of the Dirac strings. $T=2\pi/\omega_{\rm ac}$ is the time periodicity of the ac magnetic field.

To clarify the characteristics of these intrinsic collective modes, we first focus on a cubic cell with eight spins at its corners that constitute a hedgehog or an antihedgehog spin texture. The center of the (anti)hedgehog is located within this cube. In the absence of external magnetic field, the (anti)hedgehog spin texture is $C_{2x}$-, $C_{2y}$-, and $C_{2z}$-symmetric, and thus the center of mass of the solid angles spanned by four-spin plaquettes on the six cube faces coincides with the cube center. On the contrary, under application of a magnetic field along $z$-axis, the spin texture is no longer $C_{2x}$- and $C_{2y}$-symmetric, and, thereby, the center of mass of the spin solid angles deviates from the cube center along the $z$ axis.

When a magnetic field along the $z$ axis is applied to the (anti)hedgehog spin texture on a cube, the solid angles $\Omega_{\rm(A)H}^{\rm U}$ and $\Omega_{\rm(A)H}^{\rm L}$ spanned, respectively, by the four-spin plaquette on the upper face and that on the lower face of the cube become different, whereas the four-spin plaquettes on the side faces constitute solid angles of the same magnitude according to the symmetry. Consequently, the center of mass of the spin solid angles shifts along the $z$ axis as shown in Figs. 4(a) and (b). We define the displacement $\delta$ as $\delta_{\rm(A)H}/a\equiv\frac{1}{2}(|\Omega_{\rm(A)H}^{\rm U}|-|\Omega_{\rm(A)H}^{\rm L}|)/(|\Omega_{\rm(A)H}^{\rm U}|+|\Omega_{\rm(A)H}^{\rm L}|)$ with $a$ being the lattice constant.

When the magnetic field temporally oscillates along the $z$ axis, the center of mass of the spin solid angles dynamically changes its position along the $z$ axis, i.e., moves upwards ($\delta_{\rm(A)H}>0$) and downwards ($\delta_{\rm(A)H}<0$) in an oscillatory manner. The temporal displacement of the center of mass can be calculated by simulating the time evolution of localized spins when the mode is excited by ac magnetic field $\bm{h}(t)$ with the corresponding frequency. In Figs. 4(c) and 4(d), the calculated time profiles of displacements $\delta_{\rm(A)H}$ are plotted for the String B in the L2 mode and the String A in the L3 mode, respectively. We find that the hedgehogs and antihedgehogs oscillate with the same phase. Such an in-phase oscillation of the center-of-mass of each (anti)hedgehog can be regarded as coherent translational motion of the Dirac strings. Schematics of the corresponding motion are also shown in Figs. 4(e) and 4(f).

The oscillatory variation of solid angles of four-spin plaquettes orthogonal to the $z$ axis is analogous to the breathing mode in a skyrmion crystal [58]. It is known that a skyrmion crystal has one breathing mode because it is crystallographically composed of one skyrmion. In contrast, the present 4$Q$-HL states have multiple modes because it contains crystallographically inequivalent two Dirac strings with different vorticities of $\pm 1$. In this case, two localized modes, i.e., the String-A activated (L3) and the String-B activated (L2) modes are possible. The number of modes active to the longitudinal field is governed by the number of crystallographically inequivalent Dirac strings in the hedgehog-antihedgehog lattice states. This indicates that the absorption spectra in the sub-terahertz regime provide a strong clue to identifying real-space configurations of the hedgehog lattice states and their topological phase transitions caused by the monopole-antimonopole pair annihilations.

In summary, we have theoretically studied the spin-wave modes of the 4$Q$-HLs in the spin-charge coupled metallic magnets described by the chiral Kondo-lattice model. We have discovered translation modes of Dirac strings in the hedgehog lattices. It has been found that number of the collective modes is governed by the number of crystallographically inequivalent Dirac strings, which offers an experimental opportunity to reveal the spatial configurations of the hedgehogs and antihedgehogs in real magnets such as MnSi_1-xGe_x and SrFeO₃. Our findings are also expected to provide insights into dynamics of hedgehog lattices in other types of magnets including insulating systems [59, 60, 61] and their couplings to other degrees of freedom such as lattices and polarizations [62]. The nature and properties of the intrinsic excitation modes in the hedgehog lattices in magnets will open a new research field on the fundamental physics and even engineering of emergent magnetic monopoles in condensed matters.

This work was supported by JSPS KAKENHI (Grant No. 20H00337 and No. 23H04522), JST CREST (Grant No. JPMJCR20T1), and Waseda University Grant for Special Research Projects (2023C-140). R.E. was supported by a Grant-in-Aid for JSPS Fellows (Grant No. 23KJ2047). A part of the numerical simulations was carried out at the Supercomputer Center, Institute for Solid State Physics, University of Tokyo.

Supplemental Material for
“Theory of Collective Excitations in the Quadruple-Q Magnetic Hedgehog Lattices”

I Chemical potential dependence of the magnetic wavenumber

In the present study, the chemical potential is fixed at a specific value of $\mu=-3.79$ to produce a commensurate magnetic wavenumber of $Q_{\rm abs}=\pi/4$ for the quadruple-$Q$ hedgehog lattices (4$Q$-HLs) for the convenience of numerical simulations with finite-size systems. When the chemical potential varies, the shape of Fermi surafce changes, which modulates the ordering vectors $\{\bm{Q}_{\nu}\}$. In this section, we discuss that the 4$Q$-HL states emerge not only at this specific value of $\mu$ but also within a rather wide range of $\mu$ by argument based on the bare magnetic susceptibility $\chi_{\bm{q}}^{0}$. By the second-order perturbation expansion with respect to the spin-charge coupling term (the second term) of the Kondo-lattice Hamiltonian $\mathcal{H}_{\rm KLM}$ in Eq.(1) in the main text, the free energy of the system is obtained as [47],

$$F^{(2)} =-\frac{J_{\rm K}^{2}}{4}\sum_{{\bm{q}}\in{\rm BZ}}\chi_{\bm{q}}^{0}\left|{\bm{S}}_{\bm{q}}\right|^{2},$$ (S1)

with

$$\chi_{\bm{q}}^{0} =\frac{1}{N}\sum_{{\bm{k}}\in{\rm BZ}}\frac{f\left(\varepsilon_{\bm{k}}-\mu\right)-f\left(\varepsilon_{{\bm{k}}+{\bm{q}}}-\mu\right)}{\varepsilon_{{\bm{k}}+{\bm{q}}}-\varepsilon_{{\bm{k}}}},$$ (S2)

where $f(\varepsilon)$ and $\varepsilon_{\bm{k}}$ denote the Fermi distribution function and the momentum representation of kinetic energy for the first term of $\mathcal{H}_{\rm KLM}$, respectively. The maxima of $\chi_{\bm{q}}^{0}$ give the propagation vectors of magnetic helices $\{{\bm{Q}}_{\nu}\}$. Here, $\nu$ is an index of multiple equivalent peaks of $\chi_{\bm{q}}^{0}$, if any, which are connected by symmetry operations of the cubic lattice.

We know from experience that the combination of the nearest-neighbor positive electron hopping and the farther-neighbor negative hopping often induces multiple equivalent maxima of $\chi_{\bm{q}}^{0}$ on the high-symmetric lines in the Brillouin zone [47]. Here we set $t_{1}=1$ and $t_{4}=-1$, which is also presented in the main text. Fig. S1(a) shows chemical potential dependence of the ordering wavenumber at which $\chi_{\bm{q}}^{0}$ has maximal peaks. We find the wavenumber appears on a line of $q_{x}=q_{y}=q_{z}$ when $\mu$ is within a certain range of $-4.4\lesssim\mu\lesssim-2.4$. Note that the existence of the maximal peaks of $\chi_{\bm{q}}^{0}$ at such high-symmetric momenta is necessary and sufficient conditions for realization of the 4$Q$-geometry. This observation indicates that the 4$Q$-geometry is realized in the finite range of chemical potential. For further information, the Fermi-surface geometries for various values of chemical potential are presented in Fig. S1(b)-(f).

Figure S1: (a) Chemical potential dependence of the wavenumber ${\bm{q}}=(q_{x},q_{y},q_{z})$ at which the bare magnetic susceptibility $\chi_{\bm{q}}^{0}$ has maximal peaks in the crystallographic Brillouin zone. The model parameters are set as $t_{1}=1$, $t_{4}=-1$, and the system size is $N=192^{3}$. The red line represents a high-symmetric momentum path which satisfies $q_{x}=q_{y}=q_{z}$. We present only one maximum for a given chemical potential $\mu$ which satisfies $q_{x}\geq q_{y}\geq q_{z}\geq 0$ for clear visibility, even when there are multiple maxima of $\chi_{\bm{q}}^{0}$. (b)-(f) Three-dimensional views of Fermi surfaces when (b) $\mu=-4.4$, (c) $\mu=-4.0$, (d) $\mu=-3.5$, (e) $\mu=-3.0$, and (f) $\mu=-2.4$.

II Perturbational analysis for the effective fourth-order spin-spin interactions

The 4$Q$-geometry realized by the two-body spin-spin interactions discussed in the previous section is necessary but not sufficient to stabilize the 4$Q$-HL states. In order to stabilize the 4$Q$ magnetic structures, the four-body spin-spin interactions are indispensable, whose effects have been discussed in previous studies [49, 50, 47]. In this section, following a procedure presented in Ref. [47], we evaluate contributions of the effective fourth-order spin-spin interactions to the free energy by using the perturbational analysis of $\mathcal{H}_{\rm KLM}$. By sequential perturbation expansions with respect to the Kondo coupling in the weak-coupling regime of $J_{\rm K}\ll\max_{{\bm{k}}\nu}\left[\varepsilon_{{\bm{k}}\nu}-\mu\right]$ [47], the effective fourth-order spin-spin interactions are obtained as,

$$F^{(4)} =F^{(4)}_{1}+F^{(4)}_{2}+F^{(4)}_{3}+F^{(4)}_{4}+F^{(4)}_{5}+F^{(4)}_{6},$$ (S3)

$$F_{1}^{(4)} =\frac{J_{\rm K}^{4}}{16N}\sum_{\nu}(2A_{1}-A_{2})\left({\bm{S}}_{{\bm{Q}}_{\nu}}\cdot{\bm{S}}_{{\bm{Q}}_{\nu}}\right)\left({\bm{S}}_{-{\bm{Q}}_{\nu}}\cdot{\bm{S}}_{-{\bm{Q}}_{\nu}}\right),$$ (S4)

$$F_{2}^{(4)} =\frac{J_{\rm K}^{4}}{16N}\sum_{\nu}(2A_{2})\left({\bm{S}}_{{\bm{Q}}_{\nu}}\cdot{\bm{S}}_{-{\bm{Q}}_{\nu}}\right)^{2},$$ (S5)

$$F_{3}^{(4)} =\frac{J_{\rm K}^{4}}{16N}\sum_{\nu,\nu^{\prime}}4(B_{1}+B_{2}-B_{3})\left({\bm{S}}_{{\bm{Q}}_{\nu}}\cdot{\bm{S}}_{-{\bm{Q}}_{\nu}}\right)\left({\bm{S}}_{{\bm{Q}}_{\nu^{\prime}}}\cdot{\bm{S}}_{-{\bm{Q}}_{\nu^{\prime}}}\right),$$ (S6)

$$F_{4}^{(4)} =\frac{J_{\rm K}^{4}}{16N}\sum_{\nu,\nu^{\prime}}4(-B_{1}+B_{2}+B_{3})\left({\bm{S}}_{{\bm{Q}}_{\nu}}\cdot{\bm{S}}_{{\bm{Q}}_{\nu^{\prime}}}\right)\left({\bm{S}}_{-{\bm{Q}}_{\nu}}\cdot{\bm{S}}_{-{\bm{Q}}_{\nu^{\prime}}}\right),$$ (S7)

$$F_{5}^{(4)} =\frac{J_{\rm K}^{4}}{16N}\sum_{\nu,\nu^{\prime}}4(B_{1}-B_{2}+B_{3})\left({\bm{S}}_{{\bm{Q}}_{\nu}}\cdot{\bm{S}}_{-{\bm{Q}}_{\nu^{\prime}}}\right)\left({\bm{S}}_{-{\bm{Q}}_{\nu^{\prime}}}\cdot{\bm{S}}_{{\bm{Q}}_{\nu}}\right),$$ (S8)

$$F_{6}^{(4)} =\frac{J_{\rm K}^{4}}{16N}4X\left\{\left({\bm{S}}_{{\bm{Q}}_{1}}\cdot{\bm{S}}_{{\bm{Q}}_{2}}\right)\left({\bm{S}}_{{\bm{Q}}_{3}}\cdot{\bm{S}}_{{\bm{Q}}_{4}}\right)\right.$$

$$\quad\quad\quad\left.+\left({\bm{S}}_{{\bm{Q}}_{1}}\cdot{\bm{S}}_{{\bm{Q}}_{3}}\right)\left({\bm{S}}_{{\bm{Q}}_{2}}\cdot{\bm{S}}_{{\bm{Q}}_{4}}\right)+\left({\bm{S}}_{{\bm{Q}}_{1}}\cdot{\bm{S}}_{{\bm{Q}}_{4}}\right)\left({\bm{S}}_{{\bm{Q}}_{2}}\cdot{\bm{S}}_{{\bm{Q}}_{3}}\right)+{\rm c.c.}\right\},$$ (S9)

where the coefficients are given by products of noninteracting Green’s functions $G_{\bm{k}}=\left(i\omega_{p}-\varepsilon_{\bm{k}}+\mu\right)^{-1}$ as,

$$A_{1} =\frac{T}{N}\sum_{{\bm{k}},\omega_{p}}(G_{\bm{k}})^{2}G_{{\bm{k}}-{\bm{Q}}_{\nu}}G_{{\bm{k}}+{\bm{Q}}_{\nu}},$$ (S10)

$$A_{2} =\frac{T}{N}\sum_{{\bm{k}},\omega_{p}}(G_{\bm{k}})^{2}(G_{{\bm{k}}+{\bm{Q}}_{\nu}})^{2},$$ (S11)

$$B_{1} =\frac{T}{N}\sum_{{\bm{k}},\omega_{p}}(G_{\bm{k}})^{2}G_{{\bm{k}}+{\bm{Q}}_{\nu}}G_{{\bm{k}}+{\bm{Q}}_{\nu^{\prime}}},$$ (S12)

$$B_{2} =\frac{T}{N}\sum_{{\bm{k}},\omega_{p}}(G_{\bm{k}})^{2}G_{{\bm{k}}+{\bm{Q}}_{\nu}}G_{{\bm{k}}-{\bm{Q}}_{\nu^{\prime}}},$$ (S13)

$$B_{3} =\frac{T}{N}\sum_{{\bm{k}},\omega_{p}}G_{\bm{k}}G_{{\bm{k}}+{\bm{Q}}_{\nu}}G_{{\bm{k}}+{\bm{Q}}_{\nu^{\prime}}}G_{{\bm{k}}+{\bm{Q}}_{\nu}+{\bm{Q}}_{\nu^{\prime}}},$$ (S14)

$$X =\frac{T}{N}\sum_{{\bm{k}},\omega_{p}}G_{\bm{k}}G_{{\bm{k}}+{\bm{Q}}_{1}}G_{{\bm{k}}+{\bm{Q}}_{1}+{\bm{Q}}_{2}}G_{{\bm{k}}+{\bm{Q}}_{1}+{\bm{Q}}_{2}+{\bm{Q}}_{3}}.$$ (S15)

Here $\omega_{p}$ denotes the fermionic Matsubara frequency. Fig. S2 shows temperature-dependence of each coefficient given by Eqs.(S10)-(S15) and partial free energy from each scattering process given by Eqs.(S4)-(S9) for the parameter set used for the present study, i.e., $(t_{1},t_{4},\mu)=(1,-1,-3.79)$. As seen in Fig. S2(a), the coefficient $A_{2}$ has the largest absolute value than any other coefficients at low temperatures. Hence the contribution of $F_{2}^{(4)}$ to the free energy is positive and most significant at low temperatures among all the four-body spin-spin interaction terms. Combining this result about the four-body spin-spin interactions with the knowledge about the two-body spin-spin interactions discussed in the previous section, we find that an effective model of the Kondo-lattice model examined in the present study is the bilinear-biquadratic model, whose Hamiltonian is given by,

$$\mathcal{H}_{\rm eff}=2\sum_{\nu}\left[-J|{\bm{S}}_{{\bm{Q}}_{\nu}}|^{2}+\frac{K}{N}|{\bm{S}}_{{\bm{Q}}_{\nu}}|^{4}\right].$$ (S16)

The ground-state spin configuration for this effective model was studied in the previous work [32]. In this Hamiltonian, only the contributions of $F^{(4)}_{2}$ and $F^{(2)}$ at the ordering vectors $\{\bm{Q}_{\nu}\}$ ($\nu=1,2,3,4$) are taken into account in a phenomenological manner. Even with this simplification, this effective spin model reproduces the zero-field ground state of the original Kondo-lattice model, i.e., the 4$Q$-HL state.

Figure S2: Numerically evaluated temperature-dependence of (a) the coefficients in Eqs.(S10)-(S15) and (b) the partial free energies from various scattering processes of the four-body spin-spin interactions in Eqs.(S4)-(S9). Note that the magnitudes of partial free energies in (b) are normalized by $J_{\rm K}^{4}/(16N)$. Insets are magnified plots for lower-temperature region of (a) and (b), respectively. For the numerical calculations, we use a system of $N=48^{3}$ sites and take $2^{20}$ points for the Matsubara frequencies.

III Phase diagram

Figure S3: (a), (b) Detailed phase diagram of $\mathcal{H}$ in the main text as a function of the external magnetic field $H_{z}$ for (a) the chiral case with $D=0.0002$ and (b) the nonchiral case with $D=0$. The pair of numbers $(d_{\rm A},d_{\rm B})$ in each box represents the lengths of Dirac strings A and B, and $d_{\rm A}=0$ ($d_{\rm B}=0$) means the Dirac strings A (B) do not exist or the length of Dirac strings A (B) is shorter than the lattice constant. (c) Field dependence of the phase degree of freedom $\Theta$ and the ellipticity of four helices $\varepsilon_{{\bm{Q}}_{\nu}}$ constituting the 4$Q$-HL spin structure. Since the ellipticities of the four helices are equivalent for a given magnetic field, we simply denote $\varepsilon_{{\bm{Q}}_{n}}(=\varepsilon_{{\bm{Q}}_{1}}=\varepsilon_{{\bm{Q}}_{2}}=\varepsilon_{{\bm{Q}}_{3}}=\varepsilon_{{\bm{Q}}_{4}})$. The light pink, light green, light blue, and gray backgrounds show the 4$Q$-HL phase ($N_{\rm m}$=16), the 4$Q$-HL phase ($N_{\rm m}$=8), the intermediate 4$Q$ phase, and the forced ferromagnetic phase, respectively. (d) Helix of the localized spins seen along the propagation vector. The ellipticity of the helix is defined by the ratio between lengths of the short and long axes as $\varepsilon_{{\bm{Q}}_{n}}=\eta_{{\bm{Q}}_{n}}^{\rm S}/\eta_{{\bm{Q}}_{n}}^{\rm L}$.

Figures S3(a) and S3(b) show phase diagrams of $\mathcal{H}$ in the main text as a function of the external magnetic field $H_{z}$ for the chiral case with $D=0.0002$ and the nonchiral case with $D=0$, respectively. In the nonchiral case with $D=0$, the zero-field ground state is the 4$Q$-HL with $N_{\rm m}=16$, which is composed of four vortex-type Dirac strings A and four antivortex-type Dirac strings B in its magnetic unit cell. This 4$Q$-HL spin structure can be interpreted as a superposition of two helices with lefthanded chirality and two helices with righthanded chirality [35]. In other words, pairs of helices with opposite chiralities constitute the nonchiral 4$Q$-HL. This pairwise chiralities of helices result in the cancellation of net scalar chirality for the nonchiral 4$Q$-HL state even in the presence of a magnetic field. As the magnetic field increases, the Dirac strings are shortened gradually and eventually the hedgehog and antihedgehog belonging to the string A share the same unit cube and, consequently, a transition from the nonchiral 4$Q$-HL phase with $N_{\rm m}=16$ to an intermediate 4$Q$ phase occurs. It is worth mentioning that the Dirac strings A and B always have equal lengths ($d_{\rm A}=d_{\rm B}$) upon the variation of external magnetic field $H_{z}$.

In the chiral case with $D=0.0002$ in Fig. S3(a), the zero-field ground state is again the 4$Q$-HL with $N_{\rm m}=16$ composed of four strings A and four strings B in the magnetic unit cell. However, the spin structure of this 4$Q$-HL is distinct from that of the nonchiral case. The chiral 4$Q$-HL state can be regarded as a superposition of four righthanded helices and thus has a finite scalar chirality in the presence of magnetic field [32]. As the magnetic field $H_{z}$ increases, four Dirac strings A and four Dirac strings B in the magnetic unit cell are shortened gradually. Here, the shortening rate of the strings B is faster than that of the strings A. Eventually the first topological transition occurs from the 4$Q$-HL phase with $N_{\rm m}=16$ to a new 4$Q$-HL phase with $N_{\rm m}=8$ when the strings B disappear due to the hedgehog-antihedgehog pair annihilations at the strings B. Further increase of the magnetic field leads to further shortening of the strings A. Finally, the second transition occurs from the 4$Q$-HL phase with $N_{\rm m}=8$ to the intermediate 4$Q$ phase when the hedgehog and antihedgehog belonging to the string A share the same unit cube, as seen in Fig. S7(d). As further increase in the magnetic field, the intermediate 4$Q$ phase turns into the forced ferromagnetic phase [See also Fig. S7(e)] with an abrupt magnetization jump presented in Fig. S8.

The first topological transition is characterized by abrupt changes of internal degrees of freedom in the chiral 4$Q$-HL state. We focus on the phase degree of freedom $\Theta$ shown in Fig. S3(c) which is given by the sum of the translational degree of freedom of four helices $\theta_{\nu}$ $(\nu=1,2,3,4)$ as described in the main text. In the absence of magnetic field, $\Theta$ is fixed at $\pi/3$. As the magnetic field increases, $\Theta$ decreases gradually and exhibits abrupt decrease when $H_{z}$ reaches aroud the threshold field. Eventually, $\Theta$ goes to zero when $H_{z}$ slightly exceeds the threshold field. More precisely, the first topological transition occurs shortly before $\Theta$ goes to zero, which can be understood from the results of previous studies treating a spatially continuum limit [36]. We can also see that the ellipticities of the four helices $\varepsilon_{{\bm{Q}}_{n}}=\eta_{{\bm{Q}}_{n}}^{\rm S}/\eta_{{\bm{Q}}_{n}}^{\rm L}$ [Fig. S3(d)] shows an anomalous behavior near the transition point. The ellipticity usually increases as the magnetic field increases, but it decreases near the transition point.

IV Mode localizations in the hedgehog lattice phases

Figure S4: Similar plot as Fig. 3 in the main text for the hedgehog lattice phase with $N_{\rm m}=16$ at $H_{z}=0$.

Figure S5: Similar plot as Fig. 3 in the main text for the hedgehog lattice phase with $N_{\rm m}=8$ at $H_{z}=0.005$.

In this section, we present similar data as Fig. 3 in the main text for different values of external magnetic fields. Figs. S4 and S5 show the oscillation amplitudes of the modes in the $N_{\rm m}=16$ at $H_{z}=0$ and those in the $N_{\rm m}=8$ phase at $H_{z}=0.005$, respectively. In these figures, we can see almost similar spatial distributions of oscillation amplitudes as Fig. 3 in the main text for each mode, where the L1 mode is distributed throughout the magnetic unit cell and the L2 (L3) mode is localized around the strings B (A). Note that we do not have the data of the L2 mode in Fig. S5 because the L2 mode is absent in the $N_{\rm m}=8$ phase. It is also worth mentioning that the $C_{4z}$-breaking nature seen in Figs. S4(b) and S4(c) and the $C_{4z}$-preserving nature seen in Figs. S5(b) and S5(c) are closely related to the symmetry nature of spin textures discussed in Sec. VI.

V String dynamics in the L1 mode

Figure S6: (a),  (b) Time profiles of the displacements $\delta_{\rm(A)H}$ of the Bloch points in the L1 mode for (a) String A and (b) String B. (c), (d) Schematics of the oscillatory translational motion of the Dirac strings in the L1 mode for (c) String A and (d) String B. Note that $T=2\pi/\omega_{\rm ac}$ is the time periodicity of the ac magnetic field $h_{z}(t)$.

While the L1 mode is not localized around the strings as explained in the main text, we find finite oscillation amplitudes around both the Strings A and B when we activate the L1 mode by applying an ac magnetic field with the corresponding frequency. So, regardless of how meaningful the following analysis is, we can extract the dynamics of both Strings A and B in the L1 mode by using the same method for the L2 and L3 modes. Figures. S6(a) and S6(b) show calculated time profiles of displacements $\delta_{\rm(A)H}$ for the Strings A and B in the L1 mode. For both strings, $\delta_{\rm H}$ and $\delta_{\rm AH}$ show in-phase time-periodic oscillations, which indicates translational motion of the Dirac strings. Schematics of this translational motion is shown in Fig. S6(c) and S6(d). The oscillation phases of both Strings A and B in the L1 mode are opposite to those in the L2 and L3 modes. In contrast to the dynamics of Strings B (A) in the L2 (L3) mode which shifts upward (downward) when the sign of the ac magnetic field $h_{z}(t)$ is positive (negative) as argued in the main text, both the Strings A and B in the L1 mode shift downward (upward) when the sign of $h_{z}(t)$ is positive (negative).

VI Intermediate 4$Q$ phase

Figure S7: Calculated spatial configurations of the spins around the Dirac string A and spatial maps of the local scalar spin chirality on the $xy$ planes (8$\times$8-site portion of the 16$\times$16-site planes) at selected $z$ coordinates for respective phases in the phase diagram, i.e., (a) the $4Q$-HL phase with $N_{\rm m}$=16, (b) the $4Q$-HL phase with $N_{\rm m}$=8, (c) the intermediate $4Q$ (I-4$Q$) phase, and (d) the forced ferromagnetic phase.

Figure S8: Calculated magnetization profile as a function of $H_{z}$.

In the intermediate 4$Q$ phase (I-4$Q$ phase) in the phase diagram (Fig. 2 in the main text), the hedgehog and antihedgehog connected by a Dirac string A collide to be merged, and their cores share the same cubic unit cell. This phase has nonzero scalar spin chirality around the merging point although it no longer takes a quantized value. Figure S7 shows the calculated spatial spin configurations around the Dirac string A and the spatial maps of local scalar spin chirality on the $xy$ planes at selected $z$ coordinates for respective phases. We find that the nonzero scalar spin chirality indeed remains in the I-4$Q$ phase around the hedgehog-antihedgehog merging point. In Fig. 2 in the main text, we find that the L3 mode associated with the oscillation of Dirac strings A survives in the I-4$Q$ phase. This L3 mode originates from residues of the hedgehogs and antihedgehogs as a remnant of magnetic topology. The nonzero scalar spin chirality in the I-4$Q$ phase survives and takes a rather large value even on the verge of the phase boundary to the forced ferromagnetic phase with zero scalar spin chirality. This is because the phase transition between the I-4$Q$ phase and the forced ferromagnetic (FFM) phase is of strong first-order. Indeed, in Fig. S7, we can see a nearly 90-degree flop of spins from in-plane to out-of-plane directions and abrupt vanishing of scalar spin chirality when the system enters the FFM. This first-order nature of the transition can also be seen in the calculated magnetization profile as a function of $H_{z}$ in Fig.S8, which shows an abrupt jump at the transition point.

VII Remarks on the inherent $C_{4}$-symmetry breaking

In Figs. 3(b) and 3(c) in the main text, we find the breaking of $C_{4z}$ symmetry in the spatial distribution of the spin-oscillation magnitudes in the $N_{\rm m}=16$ phase. The absence of $C_{4z}$ symmetry is attributable to the breaking of $\tilde{C}_{4z}$ symmetry in the spin configuration in equilibrium despite the Hamiltonian always preserves the $C_{4z}$ symmetry. In this section, we use the term “$\tilde{C}_{4}$ symmetry” with tilde in a broader sense than we generally mean by the term “$C_{4}$ symmetry”. Specifically, while the term “$C_{4}$ symmetry” in magnetically ordered systems usually means that the spatial spin configuration is invariant under a fourfold rotation around a particular axis, “$\tilde{C}_{4}$ symmetry” in this section means merely the invariance of the absolute values of (static) spin structure factors $S_{\bm{q}}=\left(1/\sqrt{N}\right)\sum_{i}{\bm{S}}_{i}e^{i{\bm{q}}\cdot{\bm{r}_{i}}}$ under the fourfold rotation. Namely the term does not necessarily mean the invariance of the spin configuration. The magnetic ordering with broken $\tilde{C}_{4}$ symmetry appears even in the absence of magnetic field [Figs. S4(b) and S4(c)], where only one of the three $\tilde{C}_{4}$ symmetries ($\tilde{C}_{4x}$, $\tilde{C}_{4y}$, and $\tilde{C}_{4z}$) is preserved and the other two are broken. In the case of Figs. S4(a)-(c), for example, the $\tilde{C}_{4x}$ symmetry is preserved whereas the $\tilde{C}_{4y}$ and $\tilde{C}_{4z}$ symmetries are broken. Here, two axes for the broken $\tilde{C}_{4}$ symmetries are equivalent.

In the presence of external magnetic field, the situation is more complicated. This is because the spin configuration obtained by adiabatic application of magnetic field to the zero-field magnetic texture with one $\tilde{C}_{4}$-symmey and two broken $\tilde{C}_{4}$-symmetries differs depending on whether the magnetic field is applied along the $\tilde{C}_{4}$-axis for the prserved symmetry or not. In the hedgehog lattice phase with $N_{\rm m}=16$ in the presence of magnetic field, the spin configuration obtained by applying a magnetic field along the direction parallel to one of the two axes for the broken symmetries at zero field is the ground-state spin configuration because it is lower in energy than that obtained by applying a magnetic field with the same magnitude along the unique direction parallel to the axis for the preserved symmetry. Consequently, the $\tilde{C}_{4z}$ symmetry is always broken in the hedgehog-lattice phase with $N_{\rm m}=16$ in the presence of finite magnetic field along the $z$-axis because we consider only the ground-state spin configurations but not metastable ones in this work. To obtain the ground-state spin configurations, we performed an unbiased numerical analysis based on the kernel polynomial method in the present work.

In addition, the $x$-, $y$-, and $z$-axes in the absence of magnetic field are defined to be consistent with the situation in a magnetic field along the $z$-axis . Thereby, the $\tilde{C}_{4}$-symmetry axis at zero magnetic field should not be the $z$-axis but either $x$-axis or $y$-axis. Thus, the $\tilde{C}_{4z}$ symmetry is always broken in the $N_{\rm m}=16$ phase treated in the present work irrespective of the presence or absence of the magnetic field along the $z$-axis. We also note that the $\tilde{C}_{4z}$ symmetry which is broken in the $N_{\rm m}=16$ phase is abruptly recovered when the system enters the $N_{\rm m}=8$ phase [Figs. S5(b) and S5(c)]. The inherent breaking of $\tilde{C}_{4z}$ symmetry in the hedgehog-lattice spin configurations might be related with the phase degree of freedom $\Theta$. Indeed, $\Theta$ is finite in the $N_{\rm m}=16$ phase with broken $\tilde{C}_{4z}$ symmetry, while it is zero in the $N_{\rm m}=8$ phase with $\tilde{C}_{4z}$ symmetry [Fig. S3(c)]. To clarify the relation between $\Theta$ and the symmetry breaking of the hedgehog lattice spin configurations is left for future study.