Zoology of multiple-$\bm{Q}$ spin textures in a centrosymmetric tetragonal magnet with itinerant electrons

Multiple-$\bf Q$ spin order, i.e. two- or three-dimensionally modulated periodic structure characterized by a multiple number of coexisting modulation vectors $\bf Q$, has recently attracted attention as the source of rich emergent phenomenaSkXReviewFertTwo ; SkXReviewTokura ; SkXTheoryFirst . When the local magnetic moment ${\bf m}({\bf r})$ takes a noncoplanar configuration, multiple-$\bf Q$ orders are often endowed with a nontrivial topology. For example, the superposition of multiple helical modulations can be regarded as the periodic array of skyrmions (Figure 1d), which appears as a stable particle-like object with a swirling ${\bf m}({\bf r})$ texture characterized by a nonzero integer topological charge $N_{\rm sk}=\int n_{\rm sk}dxdy$MnSi ; TEMFeCoSi . Here, the topological charge density $n_{\rm sk}({\bf r})$ is defined as

$$n_{\rm sk}=\frac{1}{4\pi}{\bf n}\cdot\left(\frac{\partial{\bf n}}{\partial x}\times\frac{\partial{\bf n}}{\partial y}\right)$$

(1)

with ${\bf n}({\bf r})={\bf m}({\bf r})/|{\bf m}({\bf r})|$, and $N_{\rm sk}$ reflects how many times ${\bf n}$ wraps a unit sphere. Magnetic skyrmions are now intensively studied as a potential candidate of high-density information bitSkXReviewFertTwo or building block of natural magnonic crystal magnonic that is potentially suitable for magnetic data processing, and further explorations of nontrivial multiple-$\bf Q$ orders with exotic origin are highly demanded.

Previously, the skyrmion lattice (SkL) state with multiple-$\bf Q$ character was mostly observed in noncentrosymmetric magnets with the relativistic Dzyaloshinskii-Moriya (DM) interactionMnSi ; TEMFeCoSi ; Cu2OSeO3_Seki ; CoZnMn_First ; GaV4S8 ; Heusler . In this case, the magnetic phase diagram is rather simple, where the single-$\bf Q$ helical phase is realized in zero magnetic field and the multiple-$\bf Q$ skyrmion lattice phase is stabilized by the application of magnetic field $\bf H$. On the other hand, the skyrmion formation has recently been discovered in centrosymmetric systems, such as Gd₂PdSi₃Gd2PdSi3 , Gd₃Ru₄Al₁₂GdKagome , and GdRu₂Si₂GdRu2Si2 . Notably, these centrosymmetric Gd-compounds commonly host an extremely small skyrmion diameter ($\sim 2$ nm), which is one order of magnitude smaller than conventional DM-based noncentrosymmetric compounds. For such systems, several alternative skyrmion formation mechanisms have been proposedFrustration1 ; Frustration2 ; AkagiJPSJ ; AkagiPRL ; HayamiModel1 ; HayamiModel3 ; BatistaModel1 ; GdRu2Si2_STM ; Batista2 ; utesov , while their experimental verification has rarely been achieved.

In this work, we have exhaustively studied magnetic structures for a prototypical skyrmion-hosting centrosymmetric tetragonal magnet GdRu₂Si₂, by performing the resonant X-ray scattering (RXS) experiments. We identified a rich variety of double-$\bf Q$ magnetic structures, including the antiferroic order of meron/anti-meron-like textures with fractional local topological charge. The observed intricate magnetic phase diagram is successfully reproduced by the theoretical framework considering the four-spin interactions mediated by itinerant electrons and magnetic anisotropy. The present results will contribute to the better understanding of the novel skyrmion formation mechanism in this centrosymmetric rare-earth compound, and demonstrate that a variety of unique multiple-$\bf Q$ spin orders can be derived from a simple crystal lattice system with itinerant electrons.

Our target material GdRu₂Si₂ crystallizes into a centrosymmetric tetragonal structure with space group $I4/mmm$RRuSi2 , which consists of alternative stacking of the Gd square lattice and Ru₂Si₂ layers along the [001] axis (Figure 1a). Magnetism is governed by Gd³⁺ ion with rather isotropic magnetic moment ($S=7/2$, $L=0$). Below the magnetic ordering temperature $T_{N}$ = 46 K, the appearance of incommensurate magnetic order with the modulation vector Q = {0.22, 0, 0} has been reportedGdRu2Si2 .

In Figure 1a and 1b, the $H-T$ (temperature) magnetic phase diagrams of GdRu₂Si₂ for $\textbf{H}\parallel[001]$ and $\textbf{H}\parallel[100]$ are summarized. The phase boundaries are determined by the magnetization measurementsGdRu2Si2_PhaseDiagram , where H-induced multiple-step metamagnetic transitions are observed for both H directions (See Supplementary Note I and Supplementary Figure S2a). The phase II for $\textbf{H}\parallel[001]$ has been identified as the square skyrmion lattice state, i.e., the double-$\bf Q$ state described by the super-position of two orthogonally modulated screw structuresGdRu2Si2 and a uniform magnetization (Figure 1d). Note that the recent scanning tunneling microscopy (STM) experiment GdRu2Si2_STM revealed the two-dimensionally modulated (i.e. double-${\bf Q}$) electron density distribution in Phase II, III and possibly I, while this technique is rather sensitive to charge degree of freedom and the corresponding spin texture in the latter phases has not been confirmed directly. In addition, previous experiments mostly focused on $\textbf{H}\parallel[001]$, and the magnetic structures under $\textbf{H}\perp[001]$ remain unexplored.

To eluciadate the detailed magnetic structure in each phase, magnetic resonant X-ray scattering (RXS) experiments in resonance with the Gd L₂ edge at 5K were performed. First, we examined the development of the magnetic modulation vector ${\bf Q}$ for $\textbf{H}\parallel[100]$, by exploring the magnetic satellite reflections around a fundamental Bragg peak indexed as $(0,4,0)\pm{\bf Q}$ (Figure 2). Figure 2a and 2b indicate the line-scan profile along ($\delta$, 4, 0) and (0, 4-$\tau$, 0), which allows the identification of the magnetic modulation vectors ${\bf Q}_{1}$ = ($q$, 0, 0) and ${\bf Q}_{2}$ = (0, $q$, 0) that are parallel and perpendicular to the external magnetic field $\bf H$, respectively. Here, the sample was initially cooled at $\mu_{0}H=0$ T, and the data points represented by open symbols were obtained in a field-increasing process from 0 T to 5 T. The ones represented by closed symbols were obtained in the subsequent field-decreasing process from 5 T to 0 T. Magnetic-field dependence of the wavenumber $q$ and integrated intensity for these magnetic satellite reflections are also plotted in Figure 2c and 2d. After the initial zero-field cooling (i.e. the phase I), the magnetic satellite reflections are found for both $[100]$ and $[010]$ directions. By applying $\textbf{H}\parallel[100]$, the phase IV is stabilized for 2 T $<\mu_{0}H<$ 4 T, where the magnetic satellite reflection is observed only parallel to $\bf H$. This result indicates that the phase IV is the single-$\bf Q$ state with ${\bf Q}_{1}\parallel{\bf H}$. Above 4 T (i.e. phase III’), the magnetic satellite reflection reappears along the both directions that are parallel and perpendicular to $\bf H$, which suggests that the phase III’ is a double-$\bf Q$ state. When the $H$-value is reduced from the phase III’, only ${\bf Q}_{1}\parallel{\bf H}$ survives in the phase IV, but then both ${\bf Q}_{1}\parallel{\bf H}$ and ${\bf Q}_{2}\perp{\bf H}$ reappear in phase I. Here, the phase I is characterized by the anisotropic magnetic satellite reflections with ${\bf Q}_{1}=(0.219,0,0)$ and ${\bf Q}_{2}=(0,0.224,0)$. This indicates that the phase I is an anisotropic double-$\bf Q$ state which breaks the four-fold-symmetry. It leads to the appearance of two kinds of domains with opposite anisotropy in the phase I, and the application of $\textbf{H}\parallel[100]$ selects one of the magnetic domains. As detailed in Supplementary Note V, the phase V appearing for larger $H$-value has been identified as a single-$\bf Q$ state with ${\bf Q}_{2}\perp{\bf H}$. In case of $\textbf{H}\parallel[001]$, both ${\bf Q}_{1}$ and ${\bf Q}_{2}$ magnetic satellite peaks are always observed in the phase I, II and III as discussed in Supplementary Note VI.

Next, to understand the detailed relationship between the magnetic phases stabilized for $\textbf{H}\parallel[100]$ and $\textbf{H}\parallel[001]$, $M$-$H$ profiles have been measured for various directions of $\bf H$ rotated within the (010) plane (Supplementary Figure S2a). On the basis of these magnetization data, magnetic phase diagram as a function of $H_{[100]}$ and $H_{[001]}$ (i.e. the [100] and [001] components of $\bf H$) at 5 K is summarized in Figure 3j (See Supplementary Note II for the detail). Here, the angle between the $\bf H$-direction and the [001] axis is defined as $\theta$, and the corresponding $\theta$-dependence of magnetization at 3 T and 7 T are plotted in Figure 3k and 3l, respectively. At 3 T, the $M$-$\theta$ profile exhibits a clear step-like anomaly at $\theta=60^{\circ}$, which represents the transition between the phases III and IV (Figure 3k). In contrast, the $M-\theta$ profile at $\mu_{0}H$ = 7 T shows a smooth behavior without any anomaly (Figure 3l), demonstrating the continuous transformation between the phase III in $\textbf{H}\parallel[001]$ and the phase III’ (double-$\bf Q$ state) in $\textbf{H}\parallel[100]$. This observation suggests that the phase III is also a double-$\bf Q$ magnetic state, which is consistent with the recent STM observation of square lattice manner of charge density modulation in the phase IIIGdRu2Si2_STM . To summarize, we conclude that the phases I, II, III, and III’ are double-$\bf Q$ states, and that the phases IV and V are single-$\bf Q$ states.

To investigate the detailed spin texture in each phase, we have performed the polarization analysis of the scattered X-ray in the RXS experiments. The schematic illustration of the experimental setup is shown in Figure 3c. The magnetic field H is applied perpendicular to the scattering plane. The incident X-ray is polarized parallel to the scattering plane ($\pi$-polarization). The scattered X-ray beam may consist of two polarization components parallel ($\pi^{\prime}$) and perpendicular ($\sigma^{\prime}$) to the scattering plane, and their intensities ($I_{\pi-\pi^{\prime}}$ and $I_{\pi-\sigma^{\prime}}$) are measured separately. When the magnetic structure ${\bf m}({\bf r})$ includes the modulated spin component ${\bf m}_{\bf Q}\exp(i{\bf Q}\cdot{\bf r})+c.c.$ (with ${\bf m}_{\bf Q}$ being a complex vector), the magnetic scattering intensity $I$ is given byRXS_Rule

$$I\propto|({\bf e}_{i}\times{\bf e}_{f})\cdot{\bf m}_{\bf Q}|^{2}.$$

(2)

Here, ${\bf e}_{i}$ and ${\bf e}_{f}$ are unit vectors representing the polarization of incident and scattered X-ray beams, respectively. On the basis of the above relationship, we can separately evaluate the amplitude of each modulated spin component.

In the following, we discuss the case for the double-$\bf Q$ structure in phase III’ at $\mu_{0}H=$ 5 T applied along the [100] direction, as an example. First, the magnetic satellite peaks at $(0,4,0)\pm{\bf Q}$ are investigated with the setup shown in Figure 3a. In this configuration, the propagation vector ${\bf k}_{i}$ of incident X-ray is almost parallel to the $[0{\bar{1}}1]$ direction. Eq. (2) suggests that $I_{\pi-\pi^{\prime}}$ and $I_{\pi-\sigma^{\prime}}$ mainly reflect the [100] and $[0\bar{1}1]$ components of ${\bf m}_{\bf Q}$, respectively. Figure 3d and 3g indicate the line scan profiles of the $I_{\pi-\pi^{\prime}}$ and $I_{\pi-\sigma^{\prime}}$ intensities for ${\bf Q}_{1}$ = ($q$, 0, 0) and ${\bf Q}_{2}$ = (0, $q$, 0), respectively. In both cases, $I_{\pi-\sigma^{\prime}}$ is clearly observed, while $I_{\pi-\pi^{\prime}}$ is almost negligible. This suggests that ${\bf m}_{{\bf Q}_{1}}$ and ${\bf m}_{{\bf Q}_{2}}$ contain [010] or [001] component, but do not have the [100] component. We have further performed a similar polarization analysis for the magnetic reflections at $(0,4,4)\pm{\bf Q}$ with the setup shown in Figure 3b, where ${\bf k}_{i}$ is almost parallel to the $[0\bar{1}0]$ direction. In this case, $I_{\pi-\pi^{\prime}}$ and $I_{\pi-\sigma^{\prime}}$ mainly reflect the [100] and [010] components of ${\bf m}_{\bf Q}$, respectively. Figure 3e indicates the line scan profiles of $I_{\pi-\pi^{\prime}}$ and $I_{\pi-\sigma^{\prime}}$ for ${\bf Q}_{1}$. The observed presence of $I_{\pi-\sigma^{\prime}}$ and absence of $I_{\pi-\pi^{\prime}}$ suggest that ${\bf m}_{{\bf Q}_{1}}$ mainly consists of the [010] component. The corresponding line scan profile for ${\bf Q}_{2}$ is also shown in Figure 3h, where the absence of $I_{\pi-\pi^{\prime}}$ and $I_{\pi-\sigma^{\prime}}$ intensities demonstrates that ${\bf m}_{{\bf Q}_{2}}$ contains neither [010] nor [100] component. Since the aforementioned results for $(0,4,0)\pm{\bf Q}_{2}$ (Figure 3g) suggest that the projection of ${\bf m}_{{\bf Q}_{2}}$ along the $[0\bar{1}1]$ direction is not zero, ${\bf m}_{{\bf Q}_{2}}$ should possess the [001] component. On the basis of the above results, we conclude that the spin texture in the phase III’ for $\textbf{H}\parallel[100]$ can be approximately described as

$${\bf m}({\bf r})\propto[(0,m_{b},0)\exp(i{\bf Q}_{1}\cdot{\bf r})+(0,0,m_{c})\exp(i{\bf Q}_{2}\cdot{\bf r})+c.c.]+M_{0}(1,0,0)$$

(3)

with $M_{0}$ representing the $\bf H$-induced uniform magnetization component. The resultant anisotropic double-$\bf Q$ spin texture for phase III’ is the superposition of two orthogonally modulated sinusoidal spin components (Figure 3f and 3i) as shown in Figure 1f.

By performing similar RXS measurements with the polarization analysis, each ${\bf m}_{\bf Q}$ component has been evaluated for all the magnetic phases (i.e. phase I - V). The experimentally identified spin textures are summarized in Figure 1c-h. (See Supplementary Note III, IV, V and VI for the detail). The phase I, which has been identified as an anisotropic double-$\bf Q$ state from the aforementioned discussion for Figure 2, turns out to possess the spin texture described by

$${\bf m}({\bf r})\propto[(0,m_{b},im_{c})\exp(i{\bf Q}_{1}\cdot{\bf r})+(m_{a},0,0)\exp(i{\bf Q}_{2}\cdot{\bf r})]+c.c.,$$

(4)

i.e. the superposition of screw and sinusoidal spin modulations with ${\bf Q}_{1}=(0.219,0,0)$ and ${\bf Q}_{2}=(0,0.224,0)$, respectively (Figure 1c). On the other hand, the phase III induced by $\textbf{H}\parallel[001]$ is the double-$\bf Q$ state given by

$${\bf m}({\bf r})\propto[(0,m_{b},0)\exp(i{\bf Q}_{1}\cdot{\bf r})+(m_{a},0,0)\exp(i{\bf Q}_{2}\cdot{\bf r})+c.c.]+M_{0}(0,0,1),$$

(5)

which represents the square vortex lattice state described by the superposition of two sinusoidally modulated spin components (Figure 1e). The phases IV and V induced by $\textbf{H}\parallel[100]$ are the single-$\bf Q$ states described by

$${\bf m}({\bf r})\propto[(0,m_{b},im_{c})\exp(i{\bf Q}_{1}\cdot{\bf r})+c.c.]+M_{0}(1,0,0)$$

(6)

and

$${\bf m}({\bf r})\propto[(0,0,m_{c})\exp(i{\bf Q}_{2}\cdot{\bf r})+c.c.]+M_{0}(1,0,0),$$

(7)

representing the screw spin order with ${\bf Q}_{1}\parallel{\bf H}$ (Figure 1g) and the sinusoidally modulated fan-like spin order with ${\bf Q}_{2}\perp{\bf H}$ (Figure 1h), respectively. Notably, the experimentally identified ${\bf m}_{\bf Q}$ in the $H$-induced phases III, III’, IV and V are always normal to the $\bf H$-direction (Figure 1e-h). In general, the magnetic moments in the antiferromagnetic and helimagnetic states favor to be aligned perpendicular to the $\bf H$-direction, and the observed ${\bf m}_{\bf Q}\perp{\bf H}$ relationship is reasonable in this context. At 7 T, the rotation of magnetic field direction from the [001] axis (phase III represented by Eq. (5)) to the [100] axis (phase III’ represented by Eq. (3)) leads to the continuous rotation of ${\bf m}_{{\bf Q}_{2}}$ within the (010) plane keeping the ${\bf m}_{{\bf Q}_{2}}\perp{\bf H}$ relationship, which well explains the observed smooth $M$-$\theta$ profile in Figure 3l.

Next, we discuss the microscopic mechanism to realize such an intricate magnetic phase diagram with a rich variety of double-$\bf Q$ orders. In the following, we demonstrate that the observed phase diagram and spin textures can be well reproduced by the effective Hamiltonian given by

$$\mathcal{H}=2N\sum_{\nu=1,2}\left[-J\left(\sum_{\alpha=x,y,z}\Gamma_{{\bf Q}_{\nu}}^{\alpha}m_{{\bf Q}_{\nu}}^{\alpha}m_{-{\bf Q}_{\nu}}^{\alpha}\right)+K\left(\sum_{\alpha=x,y,z}\Gamma_{{\bf Q}_{\nu}}^{\alpha}m_{{\bf Q}_{\nu}}^{\alpha}m_{-{\bf Q}_{\nu}}^{\alpha}\right)^{2}\right]-\sum_{i}{\bf H}\cdot{\bf m}_{i},$$

(8)

which is derived from a Kondo lattice model consisting of localized magnetic moments and itinerant electrons. This theoretical framework has originally been proposed in Refs. HayamiModel1 ; HayamiModel2 , where a set of wavenumbers ${\bf Q}_{1}=(q,0,0)$ and ${\bf Q}_{2}=(0,q,0)$ are assumed to take peaks in the bare susceptibility under the tetragonal lattice symmetry. $J(\equiv 1)$ represents the amplitude of Ruderman-Kittel-Kasuya-Yosida (RKKY) interaction stabilizing the magnetic modulation with the wavevector ${\bf Q}_{\nu}$, and $K$ represents the amplitude of four-spin interaction favoring multiple-$Q$ ordersHayamiModel1 ; HayamiModel2 ; Blugel . $\Gamma_{{\bf Q}_{1}}=(\gamma_{x},\gamma_{y},\gamma_{z})$ and $\Gamma_{{\bf Q}_{2}}=(\gamma_{y},\gamma_{x},\gamma_{z})$ describe the anisotropy of associated interactions originating from the spin-orbit coupling. The last term represents the Zeeman coupling in an external magnetic field H, and ${\bf m}_{i}$ represents the local magnetic moment at the position $i$ (See Methods section for the detailed parameters and procedures for the theoretical simulation.)

Figure 4a shows the magnetization profile for $\textbf{H}\parallel[001]$ theoretically calculated by Eq. (8) with $K=0.3$. It shows a series of step-like anomalies, in good agreement with the experimental results (Supplementary Figure S2a). The corresponding $H$-dependence of individual $(m^{\alpha}_{{\bf Q}_{\nu}})^{2}$ component and the simulated spin texture at selected $H$-values are shown in Figure 4b and 4g-i, respectively. We have found that the experimentally identified spin textures for phases I (anisotropic double-$\bf Q$ state (Figure 1c)), II (double-$\bf Q$ square skyrmion lattice state (Figure 1d)) and III (double-$\bf Q$ square vortex lattice state (Figure 1e)), as well as the $\bf H$-induced transitions among these phases (i.e. I $\rightarrow$ II $\rightarrow$ III $\rightarrow$ FM (Figure 1a)), are successfully reproduced in the present simulation. The same behaviors are also confirmed for a wide range of $K$ parameter of $0.1<K\lesssim 0.5$.

Similar theoretical calculations have been also performed for $\textbf{H}\parallel[100]$, and the obtained $H$-dependence of $\bf M$ and $(m_{{\bf Q}_{\nu}}^{\alpha})^{2}$, as well as the associated spin textures at selected $H$-values, are summarized in Figure 4c-f and 4j-l. The simulation well reproduces the experimentally identified spin textures for the phases IV (single-$\bf Q$ screw state (Figure 1g)), III’ (anisotropic double-$\bf Q$ state (Figure 1f)) and V (single-$\bf Q$ fan state (Figure 1h)), while the manner of $\bf H$-induced phase transition depends on the magnitude of $K$. The present calculation predicts the successive transitions in order of I $\rightarrow$ IV $\rightarrow$ V $\rightarrow$ FM (Figure 4c) and I $\rightarrow$ III’ $\rightarrow$ V $\rightarrow$ FM (Figure 4e) for $K=0.2$ and $K=0.3$, respectively, which suggests that the free energies for phases III’ and IV are almost degenerated. In the experiment, the transition in order of I $\rightarrow$ IV $\rightarrow$ III’ $\rightarrow$ V $\rightarrow$ FM is observed (Figure 1b) in accord with theoretical prediction, implying that the fine tuning of $K$ value and/or magnetic anisotropy would allow the full reproduction of the experimental phase transition process for $\textbf{H}\parallel[100]$. The above results suggest that Eq. (8) well captures the physics behind the skyrmion formation and the related multiple-$\bf Q$ magnetism in this system.

To understand the topological aspect of each spin texture, their corresponding spatial distribution of topological charge density $n_{\rm sk}({\bf r})$ are plotted in Figure 4m-r. The double-$\bf Q$ phases I, II and III are characterized by the common texture of in-plane spin component (Figure 4g-i), which possesses two vortices and two anti-vortices in the magnetic unit cell. By further considering the out-of-plane spin component, they can be considered as a periodic lattice of meron and anti-meron-like texture, where magnetic moments in meron and anti-meron wrap a half of unit sphere and are characterized by $N_{\rm sk}=-1/2$ and $N_{\rm sk}=+1/2$, respectively Meron1 ; Meron3 ; Meron2 . According to Eq. (1), the sign of $N_{\rm sk}$ is affected by the orientation of out-of-plane spin component at their core positions and the swirling manner of in-plane spin component in the surrounding area. For example, the magnetic unit cell of the phase II consists of a core-down vortex (i.e. a meron with $N_{\rm sk}=-1/2$), a core-up vortex (i.e. an anti-meron with $N_{\rm sk}=+1/2$) and two core-up anti-vortices (i.e. two merons), which can be viewed as the square lattice of skyrmion ($N_{\rm sk}=-1$) (Figure 4n). On the other hand, the phases I and III are characterized by the different orientation of core magnetization, and can be considered as the antiferroic order of meron/anti-meron-like textures into the stripe and checkerboard patterns, respectively (Figure 4m and 4o). Likewise, the double-$\bf Q$ phase III’ can be also interpreted as the checkerboard antiferroic order of meron/anti-meron-like textures (Figure 4q). Here, the phases III and III’ smoothly transform into each other by rotating magnetic field, keeping the ${\bf m}_{\bf Q}\perp{\bf H}$ relationship. Since $n_{\rm sk}({\bf r})$ is invariant under the global rotation of magnetic moments according to Eq. (1), these two phases are characterized by the common spatial distribution of $n_{\rm sk}$ as shown in Figure 4o and 4q.

In general, conduction electrons interacting with such noncoplanar magnetic textures feel the local emergent magnetic field $b_{\rm em}^{z}({\bf r})=(h/e)n_{\rm sk}({\bf r})$ due to the additional quantum Berry phaseSkXReviewTokura ; THE ; EmergentEfield . In case of GdRu₂Si₂ with magnetic modulation period $\lambda\sim 1.9$ nmGdRu2Si2 , the amplitude of $b_{\rm em}^{z}({\bf r})$ should be in order of $\sim 100$ T. The above analysis suggests that the presently identified double-$\bf Q$ phases are characterized by the periodic modulation of giant emergent magnetic field $b_{\rm em}^{z}({\bf r})$, which will lead to the exotic manner of local electron dynamics (such as meandering motion of electrons). Note that $n_{\rm sk}({\bf r})$ is always negligible for the phase IV and V with single-$\bf Q$ character (Figure 4p,r), suggesting that the multiple-$\bf Q$ nature is essential for the appearance of local topological charge and associated emergent magnetic field.

In Figure 4s-x, we plot theoretically calculated reciprocal-space distribution of $|{\bf m}_{\bf Q}|$ for each phase. They always exhibit peak structures corresponding to $\bf Q_{1}$ and $\bf Q_{2}$. On the other hand, the higher order peak at the $\bf Q_{1}+Q_{2}$ position can be identified only in phase II (square skyrmion lattice phase), but not in the other double-$\bf Q$ phases. Since the magnetic scattering intensity in the diffraction experiment is generally proportional to $|{\bf m}_{\bf Q}|^{2}$, the absence of $\bf Q_{1}+Q_{2}$ reflection in phases I and III reported in GdRu2Si2 would be reasonable.

In this study, we have identified a variety of double-$\bf Q$ spin textures (including the antiferroic order of meron/antimeron-like textures with fractional local topological charge) for a prototype skyrmion-hosting centrosymmetric tetragonal magnet GdRu₂Si₂. The observed good agreement between the experimental and theoretical phase diagrams suggests a potential mechanism for the skyrmion formation in this system, where the magnetic interactions mediated by itinerant electrons and the associated magnetic anisotropy play a significant role. Note that several alternative theoretical frameworks based on magnetic frustration has also been proposed recently, some of which may be effectively mapped into the present model and possibly reproduce the observed phase diagrams by tuning their parameters Batista2 ; utesov . These novel theoretical frameworks predict the appearance of numerous multiple-$\bf Q$ orders consisting of various topological solitons such as merons, skyrmions and higher-order skyrmions with $N_{\rm sk}=\pm 1/2$, 1, and 2, respectively, depending on the symmetry of underlying crystal lattice and magnetic anisotropy HayamiModel1 ; HayamiModel2 ; HayamiModel3 ; BatistaModel1 . Previously, such multiple-${\bf Q}$ spin textures with fractional or higher-order topological charge have rarely been identified experimentally Meron1 ; Meron3 . The present results glimpse how rich variety of multiple-$\bf Q$ orders with unique topology/symmetry can be derived from a simple crystal lattice system with itinerant electrons, and highlight rare-earth intermetallics as a promising platform for the further search of exotic topological soliton ensembles with nontrivial functionality. The direct real-space observation of these textures and full theoretical reproduction of finite temperature phase diagram are the issue for the future study.