Role of Single-Ion Anisotropy in Stabilizing Higher-Order Skyrmion Crystals in D_3d-Symmetric Magnets

1. Introduction

Magnetic skyrmions, particle-like topological spin textures characterized by a quantized winding number [1,2,3,4], have been extensively studied owing to their fundamental importance in topological magnetism and their potential applications in spintronic devices [5,6,7,8,9,10,11,12,13]. The nontrivial topology of skyrmion spin textures gives rise to emergent electromagnetism, such as the topological Hall effect [14,15,16,17,18,19,20], and endows them with remarkable stability against perturbations [13,21,22,23,24,25]. In many noncentrosymmetric systems, skyrmion crystals (SkX), which correspond to a periodic alignment of skyrmions, are stabilized by the Dzyaloshinskii–Moriya (DM) interaction [26,27], which favors spin twisting [3,4,28,29,30,31,32,33]. This mechanism has been realized in various chiral and polar magnets, including metallic B20 compounds like MnSi [17,34,35,36,37,38], Fe_1−xCo_xSi [39,40,41,42], and FeGe [43,44], and multiferroic insulators like Cu₂OSeO₃ [45,46,47] and VOSe₂O₅ [48,49], where the interplay between ferromagnetic exchange interaction and the DM interaction leads to the formation of hexagonal and tetragonal SkXs under external magnetic fields.

In contrast, the emergence of SkXs in centrosymmetric magnets has attracted increasing attention in recent years, as it demonstrates that the DM interaction is not a prerequisite for the formation of SkX spin textures [50,51,52]. In such systems, competing symmetric exchange interactions and magnetic anisotropies can induce multiple-Q states, in which several spin density waves with distinct propagation wave vectors coexist [53,54,55,56,57]. The superposition of these spin density waves naturally gives rise to noncoplanar spin configurations with finite scalar spin chirality, thereby mimicking the topological properties of skyrmions even in the absence of DM interaction [58]. This mechanism provides an alternative and broader route toward realizing SkXs in a wide class of frustrated magnets, including both itinerant and localized spin systems on triangular and square lattices. Experimental realizations have also been reported in centrosymmetric magnets such as the rare-earth intermetallics Gd₂PdSi₃ [59,60,61,62,63], Gd₃Ru₄Al₁₂ [64,65,66,67], and GdRu₂Si₂ [68,69,70,71].

Among the various types of SkXs, higher-order SkXs, which possess a skyrmion number greater than one, have recently been recognized as intriguing topological spin textures [72,73,74,75]. These textures host multiple vortex cores within a single magnetic unit cell and exhibit intricate internal structures distinct from conventional SkX with the skyrmion number of one. The growing interest in higher-order SkXs stems from their potential to generate topological Hall and nonreciprocal transport responses even without external magnetic fields, net magnetization, or relativistic spin–orbit coupling [76,77,78]. Theoretically, such higher-order SkXs can emerge from multi-spin interactions [73,79,80,81] and bond-dependent anisotropy [82,83], whereas experimental observations of higher-order SkX in centrosymmetric materials have not yet been reported. Nevertheless, the microscopic mechanisms that determine the stability range, internal structure, and response of these higher-order SkXs remain incompletely understood. For the experimental realization of such higher-order SkXs, a further theoretical understanding of their stabilization mechanisms is crucial.

The $D_{3\mathrm{d}}$ point-group symmetry offers an ideal platform for exploring higher-order SkXs with the skyrmion number of two, as its symmetry allows characteristic anisotropic spin interactions that naturally bring about the formation of such phases [84]. In this work, we investigate a classical spin model on a two-dimensional triangular lattice possessing $D_{3\mathrm{d}}$ symmetry by focusing on the role of single-ion anisotropy on the stability of the higher-order SkX. By performing simulated annealing calculations, we construct the low-temperature zero-field phase diagram as functions of both the single-ion anisotropy and the $D_{3\mathrm{d}}$-type magnetic anisotropy. Our results reveal a broad region in parameter space where the higher-order SkX becomes energetically favorable. Moreover, our analysis reveals that the application of a magnetic field induces a series of topological phase transitions, whose nature depends sensitively on both the sign and the strength of the single-ion and $D_{3\mathrm{d}}$-type anisotropies. These findings shed light on the possibility of the higher-order SkXs in centrosymmetric $D_{3\mathrm{d}}$ magnets even without the DM interaction.

The remainder of this paper is organized as follows. In Section 2, we introduce an effective spin model on the triangular lattice with $D_{3\mathrm{d}}$ symmetry and present the relevant magnetic anisotropy terms, including the single-ion and $D_{3\mathrm{d}}$-type anisotropies. The numerical method based on simulated annealing is also described there. In Section 3, we present the ground-state phase diagram at zero magnetic field and show the stability region of the higher-order SkXs by the cooperative effects of the two anisotropies. We then discuss in detail the real-space spin textures, scalar spin chirality distributions, and spin structure factors appearing in the phase diagram. We also examine the effects of an external magnetic field and demonstrate a variety of field-induced phase transitions for representative parameter sets. Finally, Section 4 summarizes our main findings and discusses their implications for the realization of higher-order SkXs in centrosymmetric magnets.

2. Model and Method

To clarify the influence of single-ion anisotropy on the stability of higher-order SkX under centrosymmetric $D_{3\mathrm{d}}$ point-group symmetry, we consider a classical spin model with magnetic anisotropic interactions on a two-dimensional triangular lattice. The Hamiltonian is written as

$$\begin{array}{c}\hfill \mathcal{H}=-2J\sum _{\nu }\sum _{\alpha \beta }{I}_{{\mathbf{Q}}_{\nu }}^{\alpha \beta }{S}_{{\mathbf{Q}}_{\nu }}^{\alpha }{S}_{-{\mathbf{Q}}_{\nu }}^{\beta }-A\sum _i{\left(S_i^z\right)}^2-H\sum _i{S}_i^z,\end{array}$$

(1)

where ${\mathbf{S}}_i$ denotes a classical spin vector with unit length ($|{\mathbf{S}}_i|=1$) at site i. The spin components in momentum space at wave vector ${\mathbf{Q}}_{\nu }$ are defined as ${\mathbf{S}}_{{\mathbf{Q}}_{\nu }}=(1/\sqrt{N}){\sum }_i{\mathbf{S}}_i{e}^{-i{\mathbf{Q}}_{\nu }·{\mathbf{r}}_i}$, with N and ${\mathbf{r}}_i$ being the total number of spins and the position vector at site i, respectively. The exchange constant J sets the energy scale and is fixed to $J=1$ throughout this study. The second term in Eq. (1) represents the single-ion anisotropy A, originating from spin–orbit coupling under the local crystalline electric field. It controls the preferred spin orientation relative to the z axis: $A>0$ corresponds to an easy-axis anisotropy, whereas $A<0$ favors an easy-plane alignment. The single-ion anisotropy has been extensively studied in discussing the stability of the SkXs in various systems, including noncentrosymmetric systems [38,53,85,86,87,88,89]. The third term in Eq. (1) stands for the Zeeman coupling in the presence of the out-of-plane magnetic field H.

We restrict the interaction to three symmetry-equivalent ordering wave vectors,

$$\begin{array}{c}\hfill {\mathbf{Q}}_1=(\pi /3,0),{\mathbf{Q}}_2=(-\pi /6,\sqrt{3}\pi /6),{\mathbf{Q}}_3=(-\pi /6,-\sqrt{3}\pi /6),\end{array}$$

(2)

by supposing that the dominant magnetic instabilities on the triangular lattice occur at these ordering wave vectors. Such a situation arises from competing exchange interactions in frustrated insulating magnets [90,91,92,93] and from the Ruderman–Kittel–Kasuya–Yosida (RKKY) interaction [94,95,96] in itinerant metallic magnets [97] so that the instability toward the single-Q spiral state is induced [98,99,100]. All other $\mathbf{q}$ components are neglected under the assumption that the magnetic susceptibility exhibits sharp peaks at these ${\mathbf{Q}}_{\nu }$. Such a reduced description effectively captures the essential physics of multiple-Q states in a numerically efficient way and has been successfully applied to a variety of frustrated systems, including SkX and other multiple-Q states [101,102,103,104,105,106,107,108], such as Y₃Co₈Sn₄ [109], EuPtSi [110], and EuNiGe₃ [111].

The matrix $I_{{\mathbf{Q}}_{\nu }}^{\alpha \beta }$ in Eq. (1) encodes the symmetry-allowed anisotropic exchange interactions associated with the trigonal $D_{3\mathrm{d}}$ point-group symmetry. Its explicit tensor form is given by

$$\begin{array}{cc}\hfill I^{\mathrm{D}3\mathrm{d}}\equiv & I_{{\mathbf{Q}}_1}^{yz}=I_{{\mathbf{Q}}_1}^{zy}=-{\displaystyle \frac{2}{\sqrt{3}}}{I}_{{\mathbf{Q}}_2}^{xz}=-{\displaystyle \frac{2}{\sqrt{3}}}{I}_{{\mathbf{Q}}_2}^{zx}=-2I_{{\mathbf{Q}}_2}^{yz}=-2I_{{\mathbf{Q}}_2}^{zy}\hfill \\ \hfill & ={\displaystyle \frac{2}{\sqrt{3}}}{I}_{{\mathbf{Q}}_3}^{xz}={\displaystyle \frac{2}{\sqrt{3}}}{I}_{{\mathbf{Q}}_3}^{zx}=-2I_{{\mathbf{Q}}_3}^{yz}=-2I_{{\mathbf{Q}}_3}^{zy}.\hfill \end{array}$$

(3)

Here, $I^{\mathrm{D}3\mathrm{d}}$ parameterizes the $D_{3\mathrm{d}}$-type anisotropy, which favors sinusoidal spin modulations perpendicular to the propagation vectors ${\mathbf{Q}}_{\nu }$. In contrast to the compass- and Kitaev-type bond-dependent anisotropies [112,113,114,115,116,117], which favor purely in-plane sinusoidal spin modulations [82], the present $D_{3\mathrm{d}}$-type anisotropy gives rise to spin modulations whose oscillation plane involves both in-plane and out-of-plane components [58]. This anisotropy, allowed under trigonal crystal symmetry, introduces a momentum-resolved anisotropic exchange interaction that can promote noncoplanar multiple-Q spin textures like the higher-order SkX [84].

To determine the equilibrium spin configurations at low temperatures, we perform simulated annealing based on the Metropolis algorithm on $24\times 24$ lattices with periodic boundary conditions. Starting from a high temperature $T_0=1$, the system is gradually cooled according to $T_{n+1}=\tilde{\alpha }{T}_n$ with $\tilde{\alpha }$ between $0.999995$ and $0.999999$, until the final temperature reaches $T=0.0001$. Each annealing sequence consists of ${10}^5$–${10}^6$ Monte Carlo sweeps for both equilibration and measurements. Independent simulations are conducted from different random initial configurations to ensure reproducibility and to identify the ground-state phases.

The magnetic structures are characterized through both momentum-space and real-space observables. The spin structure factor for the spin component $\alpha =x,y,z$ is defined as

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

(4)

where $\mathbf{q}$ is the wave vector in the Brillouin zone. The in-plane contribution is given by $S_s^{xy}\left(\mathbf{q}\right)=S_s^x\left(\mathbf{q}\right)+S_s^y\left(\mathbf{q}\right)$. The corresponding Fourier amplitudes at ${\mathbf{Q}}_{\nu }$ are obtained from

$$\begin{array}{c}\hfill m_{{\mathbf{Q}}_{\nu }}^{\alpha }=\sqrt{{\displaystyle \frac{S_s^{\alpha }\left({\mathbf{Q}}_{\nu }\right)}{N}}}.\end{array}$$

(5)

The averaged magnetization along the field direction is given by

$$\begin{array}{c}\hfill M^z={\displaystyle \frac{1}{N}}\sum _i{S}_i^z.\end{array}$$

(6)

The topological nature of each spin configuration is evaluated via the scalar spin chirality,

$$\begin{array}{c}\hfill {\chi }^{\mathrm{sc}}={\displaystyle \frac{1}{N}}\sum _{\mu }\sum _{\mathbf{R}\in \mu }{\chi }_{\mathbf{R}}^{\mathrm{sc}},{\chi }_{\mathbf{R}}^{\mathrm{sc}}={\mathbf{S}}_i·({\mathbf{S}}_j\times {\mathbf{S}}_k),\end{array}$$

(7)

where $\mu =(u,d)$ distinguishes upward and downward triangles, and $(i,j,k)$ denote the three sites forming a triangle centered at $\mathbf{R}$ in counterclockwise order. This quantity is related to a discrete measure of the skyrmion number [118] and serves to distinguish topologically trivial and nontrivial magnetic phases in the low-temperature phase diagram.

3. Results

We show the stability of higher-order SkX by focusing on the interplay between the $D_{3\mathrm{d}}$-type anisotropic interaction $I^{\mathrm{D}3\mathrm{d}}$ and single-ion anisotropy A on the triangular lattice. First, we discuss the results at zero magnetic field in Section 3.1, where the higher-order SkX with the skyrmion number of two is relatively robust against the easy-axis single-ion anisotropy compared to the easy-plane one. Then, we discuss the field-induced phase transitions by introducing the magnetic field in Section 3.2.

3.1. Zero Magnetic Field

Figure 1 presents the low-temperature phase diagram obtained by the simulated annealing at the zero magnetic field ($H=0$), which reveals how the interplay between the single-ion anisotropy A and the $D_{3\mathrm{d}}$-type anisotropy $I^{\mathrm{D}3\mathrm{d}}$ governs the emergence of distinct magnetic phases. A prominent feature is the stabilization of the higher-order SkX with the skyrmion number of two (SkX-2) in a broad region at intermediate and large $I^{\mathrm{D}3\mathrm{d}}$ [58]. In contrast, smaller $I^{\mathrm{D}3\mathrm{d}}$ values favor a double-Q (2Q) phase, and for weak $I^{\mathrm{D}3\mathrm{d}}$ combined with strong $\left|A\right|$, the system stabilizes a single-Q (1Q) spiral phase. These results demonstrate that the $D_{3\mathrm{d}}$-type anisotropy plays a key role in stabilizing noncoplanar multiple-Q states, while the single-ion anisotropy modulates the degree of spin canting and the spatial structure in magnetic phases.

In the absence of the single-ion anisotropy ($A=0$), the SkX-2 state emerges as a robust ground state for $I^{\mathrm{D}3\mathrm{d}}\gtrsim 0.1$. The real-space spin configuration in the case of $I^{\mathrm{D}3\mathrm{d}}=0.3$ in Figure 2(a) exhibits a periodic triangular array of swirling textures, each magnetic unit cell containing two inequivalent vortex cores with different topological winding senses. The color scale represents $S_i^z$, where red and blue correspond to spins pointing up and down, respectively, while the arrows indicate the in-plane components forming a swirling vortex pattern. This configuration comprises two types of vortices: vortices centered at $S_i^z=+1$ with the vorticity $+1$, and antivortices centered at $S_i^z=-1$ with the vorticity $-2$. The asymmetry between their populations results in a net topological charge of $n_{\mathrm{sk}}=2$ per magnetic unit cell. Such a composite topological texture arises from the phase-locked sinusoidal superposition of three spin density waves with propagation ordering wave vectors ${\mathbf{Q}}_1$, ${\mathbf{Q}}_2$, and ${\mathbf{Q}}_3$, which are mutually separated by ${120}^{\circ }$ [84]. The scalar spin chirality map in Figure 3(a) shows a uniform hexagonal array of positive ${\chi }_{\mathbf{R}}^{\mathrm{sc}}$, corresponding to the vortex and antivortex centers. The overall positive average of ${\chi }_{\mathbf{R}}^{\mathrm{sc}}$ indicates a nontrivial topological state; it is noted that the state with the opposite sign of the scalar spin chirality has the same energy. The spin structure factor in Figure 4(a) displays sharp Bragg peaks at all three ${\mathbf{Q}}_{\nu }$ positions in both the in-plane and out-of-plane components, confirming that this higher-order SkX arises from the coherent superposition of triple-Q modulations. This spin configuration spontaneously shows a nonzero averaged magnetization, $M^z\ne 0$, which arises from a triple-Q sinusoidal superposition, even at zero magnetic field.

When the single-ion anisotropy A is introduced, the SkX-2 state continuously evolves without changing its topological character over a wide range of parameters. For the easy-axis case ($A>0$), as shown in Figure 2(b), the spins acquire stronger out-of-plane polarization, leading to a more pronounced z-spin modulation. According to the development of collinear tendency, the local scalar spin chirality becomes smaller, as shown in Figure 3(b). The spin structure factor correspondingly displays an increased amplitude in the out-of-plane components at the ${\mathbf{Q}}_{\nu }$ points, while the in-plane components are slightly suppressed, as shown in Figure 4(b). In contrast, for the easy-plane case ($A<0$), the out-of-plane spin component is reduced, while the in-plane spin component is enhanced, as shown by both the real-space and momentum-space plots in Figure 2(c) and Figure 4(c). Such a change has less influence on the local scalar spin chirality, as shown in Figure 3(c).

The systematic evolution of these quantities is summarized in Figure 5. As shown in Figure 5(a), the total magnetization $M_z$ tends to increase gradually with A; there is a kink structure around $A\simeq 0.2$. The scalar spin chirality ${\chi }^{\mathrm{sc}}$ in Figure 5(b) exhibits a dome-like variation, which is maximized around $A\simeq -0.2$. For the easy-axis anisotropy, ${\chi }^{\mathrm{sc}}$ is rapidly suppressed owing to the collinear tendency of the internal spin configuration. The Fourier amplitudes of in-plane and out-of-plane components in Figure 5(c) and Figure 5(d) show opposite trends with respect to A, which is understood from the role of the single-ion anisotrpy; $A>0$ tends to favor the z-directional oscillations, while $A<0$ tends to favor the $xy$-plane oscillations. These results clearly demonstrate that the degree of the easy-axis and easy-plane single-ion anisotropy continuously tunes the spin canting angle and the noncoplanarity in the SkX-2 phase while keeping the skyrmion number.

A particularly remarkable effect of the single-ion anisotropy appears in the spatial position of the skyrmion cores [119]. In the easy-plane regime ($A<0$), the skyrmion cores are located at interstitial sites of the triangular lattice, i.e., at the centers of the triangle plaquettes, as shown in Figure 2(c). With increasing A, the skyrmion cores shift toward the lattice sites, and in the easy-axis regime ($A>0$), they become centered directly on the atomic positions, as shown in Figure 2(b). This crossover from interstitial to on-site positions is discontinuous, taking place around $A\simeq 0.2$, where the kink structure appears in the magnetization shown in Figure 5(a), indicating the modification of the underlying multiple-Q interference pattern.

For weaker $D_{3\mathrm{d}}$-type anisotropy, the SkX-2 phase becomes unstable and gives way to simpler magnetic structures. At $I^{\mathrm{D}3\mathrm{d}}=0.08$ with small single-ion anisotropy, a double-Q state emerges, as indicated in Figure 2(d), Figure 3(d), and Figure 4(d). In this phase, two symmetry-related ordering wave vectors dominate the spin structure, producing a noncoplanar but spatially anisotropic pattern with density waves in terms of the scalar spin chirality. In contrast to the SkX-2 phase, this phase does not show the averaged scalar spin chirality in the whole system. Upon increasing $\left|A\right|$, single-Q spiral phases become stabilized. When the anisotropy is easy-axis, the spiral is of the vertical type, characterized by spin rotations mainly along the perpendicular plane to the $xy$ plane, whose real-space and momentum-space spin quantities are shown in Figure 2(e) and Figure 4(e), respectively. Meanwhile, for the easy-plane case, a conical spiral wave develops, in which spins mainly lies in an in-plane, although the z-spin modulation is weak but finite owing to the presence of the $D_{3\mathrm{d}}$-type anisotropy, as shown in Figure 2(f) and Figure 4(f). These results indicate that the $D_{3\mathrm{d}}$-type anisotropy is essential for stabilizing triple-Q SkX at zero magnetic field, while the single-ion anisotropy continuously tunes the balance between in-plane and out-of-plane spin modulations, controlling the noncoplanarity and the position of the skyrmion cores.

3.2. Nonzero Magnetic Fields

Next, we investigate the influence of an external magnetic field on the ground-state spin configurations discussed in the previous section and demonstrate a variety of field-induced phase transitions across four representative parameter sets.

3.2.1. Weak $D_{3\mathrm{d}}$-Type and Easy-Axis Anisotropies

We consider the case of weak $D_{3\mathrm{d}}$-type anisotropy under an easy-axis single-ion anisotropy, corresponding to $I^{\mathrm{D}3\mathrm{d}}=0.1$ and $A=0.3$, when the out-of-plane magnetic field H is applied. At zero magnetic field, the system forms the single-Q state without the net scalar spin chirality. By increasing the magnetic field, the single-Q state turns into the SkX-2 with the jumps of the magnetization [Figure 6(a)] and the scalar spin chirality [Figure 6(b)]. Although the Fourier amplitudes of the magnetic moments at ${\mathbf{Q}}_1$–${\mathbf{Q}}_3$ change in this state for $0.15\lesssim H\lesssim 0.65$, as shown in Figure 6(c) and Figure 6(d), the real-space spin configurations are almost the same, as shown in Figure 7(a)–Figure 7(d). The slight difference found in the position of the skyrmion core, which is located at the lattice (interstitial) site for small (large) H. This result indicates that the SkX-2 can appear even for small $D_{3\mathrm{d}}$-type anisotropy when the external magnetic field is applied.

As the magnetic field is further increased to $H\gtrsim 0.65$, the system undergoes a topological transition from the SkX with the skyrmion number of two to another SkX with the skyrmion number of one, followed by a topologically trivial multiple-Q state. The real-space snapshots in Figure 7(e) and Figure 7(f) show that the core position of $S_i^z<0$ changes from the vortex with the vorticity 2 to that with the vorticity 1. This results in the reduction of the skyrmion number from two to one in the magnetic unit cell. Indeed, this transition is accompanied by a sudden drop in ${\chi }^{\mathrm{sc}}$, as shown in Figure 6(b). Finally, the SkX with the skyrmion number of one changes into the topologically trivial triple-Q state without the skyrmion number. In this way, the topological transitions in terms of the skyrmion number occur as $0\to 2\to 1\to 0$ in this parameter region.

3.2.2. Weak $D_{3\mathrm{d}}$-Type and Easy-Plane Anisotropies

Next, we discuss the magnetic-field evolution for the case of weak $D_{3\mathrm{d}}$-type anisotropy under an easy-plane single-ion anisotropy, corresponding to $I^{\mathrm{D}3\mathrm{d}}=0.1$ and $A=-0.3$, where a coplanar single-Q spiral state is the ground state at zero magnetic field. Upon applying a small out-of-plane magnetic field, the magnetization is gradually developed, as shown in Figure 8(a), while the scalar spin chirality remains zero in the whole H region, as shown in Figure 8(b). This indicates that the topologically nontrivial spin configuration, i.e., SkX, does not appear in the case of the easy-plane single-ion anisotropy for weak $D_{3\mathrm{d}}$-type anisotropy.

Meanwhile, one finds that the single-Q spiral state changes into the triple-Q state close to the saturation magnetic field, as shown in Figure 8(c) and Figure 8(d). The spin configuration of this triple-Q state is characterized by the periodic alignment of vortex and antivortex without the z-spin oscillations, whose real-space spin configuration is presented in Figure 9; no net scalar spin chirality emerges in this state. A similar triple-Q state has also been found in other spin models, where the frustrated competing interactions play an important role [54,120].

3.2.3. Strong $D_{3\mathrm{d}}$-Type and Easy-Axis Anisotropies

Next, we discuss the situation with strong $D_{3\mathrm{d}}$-type anisotropy, $I^{\mathrm{D}3\mathrm{d}}=0.3$. In the case of the easy-axis single-ion anisotropy at $A=0.3$, the zero-field state corresponds to the SkX-2. The SkX-2 phase is robust against the external magnetic field up to $H\simeq 1.3$, which is in contrast to the weak $D_{3\mathrm{d}}$-type anisotropy with $I^{\mathrm{D}3\mathrm{d}}=0.1$. The spin configuration of the SkX-2 smoothly changes by increasing the magnetic field, as found by the behaviors of the magnetization in Figure 10(a) and the scalar spin chirality in Figure 10(b). For the small magnetic field around $H\simeq 0.2$, the intensities of magnetic moments at ${\mathbf{Q}}_1$–${\mathbf{Q}}_3$ are inequivalent, as shown in Figure 10(c) and Figure 10(d), although their real-space spin configurations are almost the same, as shown in Figure 11(a) and Figure 11(b), indicating the same topological property irrespective of H.

By further increasing the magnetic field, the SkX-2 turns into the SkX with the skyrmion number of one at $H\simeq 1.3$, similar to the case of $I^{\mathrm{D}3\mathrm{d}}=0.1$. The real-space spin configuration in Figure 11(c) resembles that in Figure 7(e). The SkX with the skyrmion number of one finally turns into the fully polarized state with a clear jump of all the physical quantities in Figure 10. Thus, the strong $D_{3\mathrm{d}}$-type anisotropy tends to favor the SkX-2 even under the magnetic field, and the topological transition between the SkXs with different skyrmion numbers is commonly found in the case of the easy-axis single-ion anisotropy.

3.2.4. Strong $D_{3\mathrm{d}}$-Type and Easy-Plane Anisotropies

Finally, we discuss the case of strong $D_{3\mathrm{d}}$-type anisotropy under an easy-plane single-ion anisotropy, corresponding to $I^{\mathrm{D}3\mathrm{d}}=0.3$ and $A=-0.3$, when the out-of-plane magnetic field H is applied. At zero magnetic field, the SkX-2 is stabilized, where the finite scalar spin chirality is present. As the magnetic field increases, the spins tend to align the magnetic-field direction in a gradual way while keeping the triple-Q spin configuration with equal intensities as shown in Figure 12(c) and Figure 12(d). Such a gradual change of the spin configuration is found in the continuous change of the magnetization in Figure 12(a) and the scalar spin chirality in Figure 12(b). The real-space spin configurations for different values of H are shown in Figure 13(a), Figure 13(b), and Figure 13(c), where the sign of $S_i^z$ in the skyrmion core around $S_i^z=-1$ at zero field is reversed while increasing the magnetic field. Thus, the topological phase transition occurs in the intermediate field, although the spin and scalar spin chirality quantities show a continuous change. Such a field-induced change of the spin configuration is clearly different from that for the easy-axis single-ion anisotropy, the latter of which shows a clear phase transition from the SkX-2 to the SkX with the skyrmion number of one.

4. Conclusions

We have theoretically investigated the emergence and stability of higher-order SkX in centrosymmetric magnets with the trigonal $D_{3\mathrm{d}}$ symmetry by systematically analyzing the effects of the single-ion anisotropy and the anisotropic exchange interaction that arises from the $D_{3\mathrm{d}}$ symmetry. Using real-space simulated annealing calculations of a classical spin model on the two-dimensional triangular lattice, we constructed comprehensive phase diagrams at both zero and finite magnetic fields and clarified the microscopic mechanisms responsible for stabilizing nontrivial multiple-Q topological spin textures.

At zero magnetic field, we demonstrated that the cooperative interplay between the $D_{3\mathrm{d}}$-type magnetic anisotropy and the single-ion anisotropy gives rise to a higher-order SkX with the skyrmion number of two. This SkX phase emerges as a robust ground state in a wide parameter range of both parameters. We also showed that the position of the skyrmion core undergoes a crossover from on-site to interstitial with changing the degree of single-ion anisotropy, indicating the sensitivity of the real-space spin texture to local anisotropic terms. Under finite magnetic fields, a series of field-induced topological transitions was observed. For both weak and strong $D_{3\mathrm{d}}$-type anisotropies, the system exhibits a sequential evolution of spin textures as $2\to 1\to 0$ in terms of the skyrmion number under the easy-axis single-ion anisotropy. The SkX with the skyrmion number of two remains remarkably stable over a wide field range, particularly in the presence of strong $D_{3\mathrm{d}}$-type anisotropy. In contrast, under easy-plane anisotropy, the higher-order SkX only appears for the strong $D_{3\mathrm{d}}$-type anisotropy, which turns into the topologically trivial triple-Q state without the skyrmion number. These results indicate that the $D_{3\mathrm{d}}$-type anisotropy is essential for the formation and stability of higher-order SkXs, while the single-ion anisotropy affects the topological phase transitions driven by the magnetic field. Our findings provide microscopic insight into the stabilization mechanism of higher-order SkXs in centrosymmetric trigonal magnets, suggesting that the interplay between local and momentum-resolved anisotropies can realize complex noncoplanar spin textures even in the absence of the DM interaction.