Quadruple-Q Skyrmion Crystal in Centrosymmetric Body-Centered Tetragonal Magnets

Abstract

We conduct a numerical investigation into the stability of a quadruple-Q skyrmion crystal, a structure generated by the superposition of four spin density waves traveling in distinct directions within three-dimensional space, hosted on a centrosymmetric body-centered tetragonal lattice. Using simulated annealing applied to an effective spin model that includes momentum-resolved bilinear and biquadratic interactions, we construct a magnetic phase diagram spanning a broad range of model parameters. Our study finds that a quadruple-Q skyrmion crystal does not emerge within the phase diagram when varying the biquadratic interaction and external magnetic field. Instead, three distinct quadruple-Q states with topologically trivial spin textures are stabilized. However, we demonstrate that the quadruple-Q skyrmion crystal can become the ground state when an additional high-harmonic wave–vector interaction is considered. Depending on the magnitude of this interaction, we obtain two types of quadruple-Q skyrmion crystals exhibiting the skyrmion numbers of one and two. These findings highlight the emergence of diverse three-dimensional multiple-Q spin states in centrosymmetric body-centered tetragonal magnets.

1. Introduction

The skyrmion was initially introduced by Tony Skyrme in particle physics as a topologically stable object, which is characterized by a topological integer that remains invariant under continuous transformations [1,2]. In condensed matter physics, such topologically nontrivial objects have been explored in various contexts, such as crystalline liquids [3,4,5,6,7], quantum Hall magnets [8,9,10,11,12,13,14], and Bose–Einstein condensates [15,16,17,18,19,20]. Among these systems, magnetic skyrmion crystals (SkXs) featuring swirling spin textures in helical magnets have attracted significant attention in both theoretical and experimental studies [21,22,23,24,25]. The SkXs arise as superpositions of multiple spin density waves and have been identified in diverse materials. In noncentrosymmetric lattice structures, the Dzyaloshinskii–Moriya (DM) interaction plays a key role in their stabilization [23,26,27,28,29,30,31]. In centrosymmetric lattice structures, frustrated exchange interactions and/or multi-spin interactions are crucial [32,33,34]. Representative materials in the former category include MnSi [35,36], ${\mathrm{Fe}}_{1-x}{\mathrm{Co}}_x$Si [37,38], ${\mathrm{Cu}}_2{\mathrm{OSeO}}_3$ [39,40], ${\mathrm{Co}}_8{\mathrm{Zn}}_8{\mathrm{Mn}}_4$ [41,42], and EuPtSi [43,44]. Examples in the latter category include ${\mathrm{Gd}}_2{\mathrm{PdSi}}_3$ [45,46], ${\mathrm{Gd}}_3{\mathrm{Ru}}_4{\mathrm{Al}}_{12}$ [47,48], ${\mathrm{GdRu}}_2{\mathrm{Si}}_2$ [49,50], NiMnIn [51], and NiMnGa [51]. Materials hosting SkXs are promising candidates for realizing giant emergent electromagnetic responses like the topological Hall and Nernst effects [36,52,53,54,55,56,57,58,59], as well as robust current-induced motions [60,61,62,63,64,65,66,67].

SkXs can be classified into several types based on their internal degrees of freedom. For instance, the helicity and vorticity degrees of freedom give rise to Néel-type, Bloch-type, and anti-type SkXs [24]. These variations are closely associated with the DM interaction, whose form depends on the types of noncentrosymmetric lattice structures [21,22,68,69]. Another distinguishing feature of SkXs is the phase degree of freedom [70,71]. The relative phases among the constituent spin density waves influence spin configurations, resulting in structures such as the meron–antimeron crystals [45,72,73,74,75], tetra-axial vortex crystals [70], and SkXs with high skyrmion numbers [76,77], thereby altering their topological properties. Furthermore, the number of constituent waves also influences the SkX structure. For example, a double-Q superposition of the spiral waves produces a square SkX, while a triple-Q superposition results in a triangular SkX. These constituent spin density waves are often determined by the rotational symmetry of the two-dimensional lattice structures. Specifically, the double-Q wave vectors in the square SkX correspond to the fourfold rotational symmetry of the square lattice, while the triple-Q wave vectors in the triangular SkX correspond to the threefold rotational symmetry of the triangular lattice. Recent studies have further demonstrated the possibility of multiple-Q SkXs in three-dimensional systems, such as the sextuple-Q SkX on the cubic lattice [78].

In this study, we investigate the emergence of multiple-Q SkXs driven by three-dimensional ordering wave vectors. We focus on a centrosymmetric body-centered tetragonal lattice and examine an effective spin model using simulated annealing techniques. The choice of the body-centered tetragonal lattice is owed to its simple lattice structure, where the DM and frustrated interactions play a lesser role in stabilizing the SkX, although materials hosting SkXs have been found, such as ${\mathrm{GdRu}}_2{\mathrm{Si}}_2$ [49,50,79,80,81]. Our results reveal a variety of instabilities, including in SkXs and other multiple-Q states. We identify three types of quadruple-Q states that arise from distinct superpositions of ordering wave vectors, depending on the strength of the four-spin biquadratic interaction and the external magnetic field. We further demonstrate the appearance of a quadruple-Q SkX with a skyrmion number of one when high-harmonic wave–vector interactions are introduced. Additionally, a quadruple-Q SkX with a skyrmion number of two emerges when these interactions become comparable in strength to the dominant exchange interaction. These findings suggest that magnetic systems with three-dimensional ordering wave vectors provide a versatile platform for realizing diverse and topologically nontrivial spin textures.

The structure of this paper is as follows: In Section 2, we introduce the effective spin model defined on a body-centered tetragonal lattice and outline the numerical approach based on simulated annealing. In Section 3, we present the magnetic phase diagram obtained by varying the biquadratic interaction and external magnetic field. We demonstrate the emergence of three distinct quadruple-Q states, although none exhibit skyrmion characteristics. We then show that the introduction of high-harmonic wave–vector interactions stabilizes two quadruple-Q SkXs with skyrmion numbers of one and two. Finally, Section 4 provides a summary of our findings and conclusions.

2. Model and Method

We focus on a three-dimensional body-centered tetragonal lattice characterized by the primitive translational vectors ${\mathit{a}}_1=(a,0,0)$, ${\mathit{a}}_2=(0,a,0)$, and ${\mathit{a}}_3=(a/2,a/2,c/2)$, where we set $a=c=1$ without loss of generality. The space group of this lattice is $I4/mmm$. The lattice structure is shown in the right panel of Figure 1. The body-centered tetragonal lattice serves as a representative structure for realizing SkX phases, as observed in materials like ${\mathrm{GdRu}}_2{\mathrm{Si}}_2$ [49,50]. Based on this framework, the effective classical spin model is expressed as

$$\begin{array}{cc}\hfill \mathcal{H}=& -2J\sum _{\nu }({\alpha }_{\nu }{\mathit{S}}_{{\mathit{Q}}_{\nu }}·{\mathit{S}}_{-{\mathit{Q}}_{\nu }})+{\displaystyle \frac{2K}{N}}\sum _{\nu }{({\alpha }_{\nu }{\mathit{S}}_{{\mathit{Q}}_{\nu }}·{\mathit{S}}_{-{\mathit{Q}}_{\nu }})}^2-H\sum _i{S}_i^z,\hfill \end{array}$$

(1)

where ${\mathit{S}}_{{\mathit{Q}}_{\nu }}=(S_{{\mathit{Q}}_{\nu }}^x,S_{{\mathit{Q}}_{\nu }}^y,S_{{\mathit{Q}}_{\nu }}^z)$ and ${\mathit{S}}_i=(S_i^x,S_i^y,S_i^z)$ represent the localized spin vectors in momentum and real spaces, respectively, with the two being related via a Fourier transformation. The subscripts i in ${\mathit{S}}_i$ and ${\mathit{Q}}_{\nu }$ in ${\mathit{S}}_{{\mathit{Q}}_{\nu }}$ stand for the site and wave vector indices, respectively, both ranging from 1 to N, where N represents the number of sites. We set the spin length of each vector ${\mathit{S}}_i$ to unity, i.e., $|{\mathit{S}}_i|=1$.

The first term represents the bilinear exchange interaction in momentum space, whose real-space counterpart is given by pairs of two-spin interactions ${\mathit{S}}_i·{\mathit{S}}_j$. We set $J=1$ as the energy unit of the model and introduce ${\alpha }_{\nu }$ to account for the variations in interactions across different wave–vector channels. We suppose that the interactions at a few particular ordering wave vectors play a dominant role in determining the low-temperature phase diagram close to the ground state. This approximation is valid when Fermi surface nesting occurs at the relevant wave vectors, as in the case of the Ruderman–Kittel–Kasuya–Yosida (RKKY) interaction within the itinerant electron model [82,83,84]. We choose the three-dimensional ordering wave vectors of ${\mathit{Q}}_1=(\pi /5,0,\pi /5)$, ${\mathit{Q}}_2=(0,\pi /5,\pi /5)$, ${\mathit{Q}}_3=(\pi /5,0,-\pi /5)$, and ${\mathit{Q}}_4=(0,\pi /5,-\pi /5)$ by setting ${\alpha }_{\nu }=1$ for $\nu =1$–4. Additionally, we include the effects of interactions at higher-harmonic wave vectors of ${\mathit{Q}}_1$–${\mathit{Q}}_4$, since it tends to lead to the multiple-Q instability [85]; we set ${\mathit{Q}}_5={\mathit{Q}}_1+{\mathit{Q}}_4={\mathit{Q}}_2+{\mathit{Q}}_3$ and ${\mathit{Q}}_6=-{\mathit{Q}}_1+{\mathit{Q}}_2=-{\mathit{Q}}_3+{\mathit{Q}}_4$ with ${\alpha }_5={\alpha }_6=\alpha <1$. The prefactor of 2 denotes the contributions from $-{\mathit{Q}}_1$ to $-{\mathit{Q}}_6$. Interactions at other wave vectors within the Brillouin zone are neglected, as their influence on the ground-state energy is minimal when the magnetic states are predominantly governed by the superposition of ${\mathit{Q}}_1$–${\mathit{Q}}_4$. For simplicity, we also disregard the effect of symmetry-allowed magnetic anisotropies under the $I4/mmm$ symmetry, such as the Ising-type and bond-dependent anisotropies, although such an effect often gives rise to instability toward the SkX, as found in the case of the easy-axis Ising-type anisotorpy [86,87]. It is noted that the DM interaction, which tends to favor the SkX [28,88], is absent due to the centrosymmetric nature of the lattice.

The second term represents the biquadratic exchange interaction with the coupling constant K, whose real-space counterpart is generally given by four-spin interactions $({\mathit{S}}_i·{\mathit{S}}_j)({\mathit{S}}_k·{\mathit{S}}_l)$. This interaction originates from the higher-order RKKY interaction in itinerant electron systems, which can be derived through perturbation analysis with respect to the spin-charge coupling [89,90]. Since this interaction increases the energy of the single-Q spiral state compared to the multiple-Q state, it serves as a driving force for multiple-Q states, such as SkX, vortex, and bubble states [91], in centrosymmetric itinerant magnets. These states have been observed in materials like ${\mathrm{UCu}}_5$ [92], ${\mathrm{Y}}_3{\mathrm{Co}}_8{\mathrm{Sn}}_4$ [93], and ${\mathrm{GdRu}}_2{\mathrm{Si}}_2$ [50]. The final term in Equation (1) stands for the Zeeman term with the magnitude H, which accounts for the effect of an external magnetic field applied along the z-direction.

We employ simulated annealing to search for the ground-state spin configuration of the model described by Equation (1). For each set of parameters $(K,H,\alpha )$, the following procedure is carried out independently: First, we initialize the system at high temperatures $T_0=1$–10, with a random initial spin configuration. The temperature is then gradually decreased according to $T_{n+1}=0.999999T_n$ after each Monte Carlo sweep, where local spins are updated one by one using the Metropolis algorithm, where the transition probability is determined by the temperature and the change in the energy during the update. Here, $T_n$ represents the temperature at the nth-step. The final temperature is set to $T=0.001$, which is much smaller than the energy scale of the model with $J=1$. Once the system reaches the final temperature, we perform between ${10}^5$ and ${10}^6$ Monte Carlo sweeps to compute physical quantities, as outlined below. The system size is chosen to be $N={10}^3$, and periodic boundary conditions are applied. We confirmed that the following results are not qualitatively altered for larger system sizes, such as $N={20}^3$.

The magnetic phases obtained at the final temperature are distinguished by the spin structure factor and scalar spin chirality. The spin structure factor for $\eta =x,y,z$ is given by

$$\begin{array}{cc}\hfill S_s^{\eta \eta }\left(\mathit{q}\right)& ={\displaystyle \frac{1}{N}}\sum _{i,j}{S}_i^{\eta }{S}_j^{\eta }{e}^{i\mathit{q}·({\mathit{r}}_i-{\mathit{r}}_j)},\hfill \end{array}$$

(2)

where ${\mathit{r}}_i$ is the position vector at site i and $\mathit{q}$ is the wave vector in the first Brillouin zone. The uniform magnetization along the field direction is given by $M=(1/N){\sum }_i{S}_i^z$.

Meanwhile, the scalar spin chirality is defined as

$$\begin{array}{cc}\hfill {\chi }^{\mathrm{sc}}& ={\displaystyle \frac{1}{N}}\sum _i\sum _{\delta ,{\delta }^{\prime }=\pm 1}\delta {\delta }^{\prime }{\mathit{S}}_i·({\mathit{S}}_{i+\delta {\widehat{d}}_1}\times {\mathit{S}}_{i+{\delta }^{\prime }{\widehat{d}}_2}),\hfill \end{array}$$

(3)

where ${\widehat{d}}_1$ (${\widehat{d}}_2$) represents a shift by a in the $\left[100\right]$ ($\left[010\right]$) direction. The SkX phases are characterized by a nonzero value of ${\chi }^{\mathrm{sc}}$. We also calculate the scalar spin chirality at site i, which is defined as ${\chi }^{\mathrm{sc}}=(1/N){\sum }_i{\chi }_i^{\mathrm{sc}}$.

3. Results

3.1. Without the High-Harmonic Wave–Vector Interaction

Initially, we neglect the influence of high-harmonic wave–vector interactions, setting $\alpha =0$. We present the ground-state phase diagram in Figure 1, which reveals five distinct magnetic phases, excluding the fully polarized state that emerges for $H\ge 2$. None of these phases exhibit a uniform scalar spin chirality, implying that the SkX phase is absent when $\alpha =0$.

When the biquadratic interaction K is set to zero, the single-Q spiral state, which is denoted as $1Q$ in the phase diagram, appears irrespective of H. In this state, the spins exhibit a conical spiral arrangement confined to the $xy$-plane. The ordering wave vector can be selected from any of the four vectors, ${\mathit{Q}}_1$–${\mathit{Q}}_4$. The real-space spin configuration of the 1Q state is illustrated in Figure 2a. As we introduce a nonzero biquadratic interaction K, the system undergoes a transformation in the low-field regime, where the 1Q state evolves into a more complex double-Q state, denoted as $2Q^{\prime }$. The new state is characterized by a mixture of two distinct modulation patterns. One component remains as a conical spiral at one of the ordering wave vectors, while the second is a sinusoidal variation of the z-spin component at another ordering wave vector, with the first modulation being significantly stronger than the second. This dual modulation structure is what distinguishes the $2Q^{\prime }$ state. The prime symbol in $Q^{\prime }$ emphasizes that this state exhibits anisotropic multiple-Q features, where the intensities of the spin structure factors from different ordering wave vectors are not equal. For this state, the two chosen wave vectors are picked from the set ${\mathit{Q}}_1$–${\mathit{Q}}_4$. We display the real-space spin configuration of the 2$Q^{\prime }$ state in the case of the dominant conical spiral modulation at ${\mathit{Q}}_3$ and the sub-dominant sinusoidal modulation at ${\mathit{Q}}_4$ in Figure 2b. Since the spin configuration of the 2$Q^{\prime }$ state is continuously connected to that of the 1Q state by decreasing the intensity of the sub-dominant sinusoidal modulation, the phase transition from the 2$Q^{\prime }$ state to the 1Q state against H is characterized by the second-order phase transition, as shown by the magnetization process in Figure 3a.

As the biquadratic interaction K is increased, a new magnetic phase, termed the 4$Q^{\prime }$ I state, emerges in the intermediate-field region, as shown in Figure 1. The spin configuration of this phase is depicted in Figure 2c, where the bubble structure characterized by $S_i^z=-1$ appears in each $xy$ plane. Since there is no vorticity around the bubble structure in terms of the $xy$-spin component, this differs from the skyrmion; there is no emergence of uniform scalar spin chirality in this phase. The transition between the 2$Q^{\prime }$ and 4$Q^{\prime }$ I states is found to be a first-order phase transition, as evidenced by the abrupt change in the spin structure factor. On the other hand, the transition between the 1Q state and the 4$Q^{\prime }$ I state is smoother, corresponding to a second-order phase transition, as indicated by the continuous change in the magnetization curve shown in Figure 3a.

In the low-field region, the 2$Q^{\prime }$ state is replaced by the 4Q I state when the biquadratic interaction K exceeds a critical value of approximately $K\ge 0.13$. The defining feature of the 4Q I state is the equal intensity of the spin structure factor at the four different ordering wave vectors at ${\mathit{Q}}_1$–${\mathit{Q}}_4$. The real-space spin configuration of the 4Q I state, shown in Figure 4a, is nontrivial. It consists of vortex–antivortex pairs that appear in each of the $xy$ planes. For instance, in the $z=0$ plane, two distinct vortex and antivortex pairs are present, giving rise to a characteristic local magnetic structure. These vortex–antivortex pairs for different $z\ne 0$ planes are related to the ordering wave vectors. The intricate arrangement of these vortices in the $xy$ planes leads to noncoplanar spin textures, which are essential for the formation of local scalar spin chirality, as shown in Figure 5a. It is noted that there is no uniform scalar spin chirality in the whole system. The phase transitions from the 4Q I state to the 4$Q^{\prime }$ I and $1Q$ states are of first-order, the former of which are shown in Figure 3b.

Meanwhile, in the high-field region, the system undergoes another significant transformation. The 1Q state is replaced by another quadruple-Q state denoted as the 4Q II state, as shown in Figure 1. As shown by the real-space spin configuration in Figure 4b, the z-spin oscillation almost vanishes, which indicates that the spin configuration is characterized by a superposition of quadruple-Q fan states with the $xy$-spin oscillations at ${\mathit{Q}}_1$–${\mathit{Q}}_4$. The scalar spin chirality aligns in a checkerboard manner, which reflects the periodic alignment of the vortex and antivortex in terms of the $xy$-spin component, as shown in Figure 5b [94]. The phase transition between the 4Q I and 4Q II states is of first-order with the jump in the magnetization, as shown in Figure 3c.

3.2. With the High-Harmonic Wave–Vector Interaction

We next explore the influence of the high-harmonic wave–vector interaction, which is introduced as a mechanism to stabilize the SkX. The resulting magnetic phase diagram as a function of $\alpha$ and H is shown in Figure 6. We here choose $K=0.2$, since the instabilities toward the quadruple-Q states (4Q I and 4Q II states) are found for $K\ge 0.2$. As $\alpha$ varies, four additional phases emerge, including two distinct SkX phases. For $\alpha \gtrsim 0.1$, the SkX I phase develops between the 4Q I and 4Q II phases. In contrast, the 4$Q^{\prime }$ II phase is stabilized in the low-field region for $\alpha \gtrsim 0.28$, preceding the appearance of the SkX I phase. At larger values of $\alpha$, the 2Q and SkX II phases appear for $\alpha \gtrsim 0.7$ and $\alpha \gtrsim 0.9$, respectively. The detailed characteristics of the spin and scalar chirality configurations for each phase are discussed below.

For the 4$Q^{\prime }$ II state, the spin configuration in real space is illustrated in Figure 7a. Similarly to the 4$Q^{\prime }$ I state, the bubble appears in each $xy$ plane. However, unlike skyrmions, it does not carry out nonzero winding numbers. In other words, the 4$Q^{\prime }$ II state is a topologically trivial state without a net scalar spin chirality. We also show the real-space distribution of the scalar spin chirality in Figure 8a. In terms of the spin structure factor, the four intensities at ${\mathit{Q}}_1$–${\mathit{Q}}_4$ are split into two distinct intensities. The magnetic phase transition between the 4Q I and 4$Q^{\prime }$ II states is identified as first-order, accompanied by discontinuous changes in the intensities of the spin structure factor, although the magnetization varies continuously, as illustrated in the case of $\alpha =0.4$ in Figure 9a.

With the increase in H in the 4$Q^{\prime }$ II state, SkX I appears in the phase diagram in Figure 6. The spin configuration is characterized by the equal intensity of the spin structure factor at ${\mathit{Q}}_1$–${\mathit{Q}}_4$. The real-space spin configuration is given in Figure 7b, where the skyrmion core at $S_i^z=-1$ is found around $z=4$. Owing to the positive winding number around the skyrmion core, the scalar spin chirality around the skyrmion core becomes negative, as shown in Figure 8b, which results in the skyrmion number of $-1$. It is noted that the antiskyrmion crystal with a skyrmion number of $+1$ has the same energy as the SkX with a skyrmion number of $-1$, since there is no magnetic anisotropy in the present model [95]. The region of the SkX I phase extends as $\alpha$ increases, which indicates that the high-harmonic wave–vector interaction plays a significant role in stabilizing the SkX spin texture. The phase transitions between the SkX I and other magnetic phases are characterized by the first-order transitions owing to the different skyrmion numbers. The magnetization process shows a jump structure below and above the SkX I phase, as shown in Figure 9a.

We find another SkX phase (SkX II phase) when the high-harmonic wave–vector interaction is almost comparable to the dominant interaction, i.e., $\alpha \sim 1$. This spin configuration in real space is shown in Figure 10a, which includes two vortices with the positive winding number ($+1$) around the region with the negative z-spin polarization and two antivortices with the negative winding number ($-1$) around the region with the positive z-spin polarization in each $xy$ plane. In other words, the SkX II phase consists of two vortices with the skyrmion number $-1/2$ and two antivortices with the skyrmion number $-1/2$, which results in the skyrmion number of two in the magnetic unit cell. Such a feature is found in the real-space scalar spin chirality in Figure 11b. According to the high skyrmion number, the scalar spin chirality in the SkX II is larger than that in the SkX I, as shown in Figure 9b. This square SkX II phase has also been found in different spin models with bond-dependent magnetic anisotropy [96,97]. The SkX II turns into the SkX I by applying the magnetic field with a jump of the skyrmion number (scalar spin chirality) as well as the magnetization, as shown in Figure 9b.

In the high-field regime, the 2Q state emerges as the stable phase, replacing the 4Q II state. This phase is predominantly characterized by the superposition of two in-plane spin density waves with identical intensities. Figure 10b illustrates the real-space spin structure of this state, where the z-component of the spin exhibits oscillations due to contributions from high-harmonic wave vectors. The noncoplanar spin arrangement leads to the formation of a scalar spin chirality density wave, with its periodicity determined by the high-harmonic wave vectors, as depicted in Figure 11b. As the magnetic field strength increases, this state gradually changens into the fully polarized phase.

4. Summary

We have explored the potential emergence of SkXs formed by three-dimensional ordering wave vectors in the centrosymmetric body-centered tetragonal-lattice system. By conducting numerical simulations using a simulated annealing approach for an effective spin model incorporating bilinear and biquadratic interactions, we have identified key factors for stabilizing SkXs and other multiple-Q magnetic states. Initially, we demonstrated that three distinct quadruple-Q states appear in the phase diagram as the biquadratic interaction and magnetic field are varied, although no SkX phase emerges under these conditions. However, we found that introducing high-harmonic wave–vector interactions leads to the stabilization of two types of SkXs with skyrmion numbers of one and two, depending on the interaction strength. These findings suggest that three-dimensional ordering wave vectors can generate a rich variety of topological spin textures through the superposition of multiple spin density waves propagating along different spatial directions.