Skyrmion crystals in centrosymmetric itinerant magnets without horizontal mirror plane.

We theoretically investigate a new stabilization mechanism of a skyrmion crystal (SkX) in centrosymmetric itinerant magnets with magnetic anisotropy. By considering a trigonal crystal system without the horizontal mirror plane, we derive an effective spin model with an anisotropic Ruderman-Kittel-Kasuya-Yosida (RKKY) interaction for a multi-band periodic Anderson model. We find that the anisotropic RKKY interaction gives rise to two distinct SkXs with different skyrmion numbers of one and two depending on a magnetic field. We also clarify that a phase arising from the multiple-Q spin density waves becomes a control parameter for a field-induced topological phase transition between the SkXs. The mechanism will be useful not only for understanding the SkXs, such as that in Gd[Formula: see text]PdSi[Formula: see text], but also for exploring further skyrmion-hosting materials in trigonal itinerant magnets.

Introduction

A magnetic skyrmion, which is characterized by a topologically nontrivial spin texture[1-3], has been extensively studied in condensed matter physics since the discovery of the skyrmion crystal (SkX) in chiral magnets[4-6]. The SkX exhibits a nonzero topological winding number called the skyrmion number , which is defined as , where is a skyrmion density related to the solid angle consisting of three spins , , and on the triangle R: [7]. The study of the SkX has attracted much attention, as the swirling topological magnetic texture owing to nonzero gives rise to an emergent electromagnetic field through the spin Berry phase and results in intriguing transport phenomena and dynamics[8-12], such as the topological Hall effect[13,14] and the skyrmion Hall effect[15,16].

The SkXs are expressed as a superposition of three spin density waves (triple-Q state) aswhere and are the unit vectors along the in-plane and z directions, respectively. , and where and are phases of each spin density wave. A variety of the SkXs are described by Eq. (1); a superposition of spiral waves for () and or represents the Bloch-type (Néel-type) SkX, while that for and represents the anti-type SkX. The real-space spin texture for the Bloch-type SkX is shown in Fig. 1a. All the SkXs have the skyrmion number of one, , in the magnetic unit cell and breaks the spatial inversion symmetry irrespective of and [8]. We call them the SkXs. The SkXs are stabilized by the Dzyaloshinskii-Moriya (DM) interaction[17,18] in chiral/polar magnets[4,19] or the competing exchange interactions in frustrated magnets[20-22].

Figure 1

Schematic pictures of the SkXs and the crystal structure. (a–c) SkXs characterized by three spiral and sinusoidal waves along the , , and directions: (a) the SkX for , (b) the SkX for , and (c) the T-SkX for in Eq. (1). (d) Centrosymmetric trigonal structure without the horizontal mirror plane. The blue spheres represent magnetic sites, while the gray spheres shifted by () from the center of the downward (upward) triangles on the magnetic layer represent nonmagnetic sites on a layer A (B).

Meanwhile, the spiral density waves are not necessarily for the formation of the SkX. By considering the superposition of the sinusoidal waves characterized by a different , another type of the SkX can emerge, as shown in Fig. 1b[23,24]. In contrast to the SkX, this spin texture exhibits the skyrmion number of two in a magnetic unit cell ( SkX), whose spatial inversion and/or sixfold rotational symmetries are broken depending on on a discrete lattice. For example, the SkX with shown in Fig. 1b has the inversion symmetry, but the SkX with , , and shows the inversion symmetry breaking. Although the SkX seems to be rare compared to the one, it is stabilized by a multi-spin interaction in itinerant magnets [23,24] or an anisotropic symmetric exchange interaction in frustrated magnets[25]. Moreover, an isolated skyrmion with is nucleated in frustrated magnets[26,27].

In the present study, we report our theoretical discovery of the SkXs by focusing on a magnetic anisotropy that arises from the absence of the mirror symmetry in the crystal structure. By constructing a microscopic effective spin model and performing simulated annealing for triangular itinerant magnets, we show that an anisotropic Ruderman–Kittel–Kasuya–Yosida (RKKY) interaction [28-30] arising from the absence of the mirror symmetry on a magnetic layer [Fig. 1d] induces the SkXs with and . The anisotropic RKKY interaction stabilizes the SkXs even without the DM, competing exchange, and multi-spin interactions[25,31-35]. The obtained SkXs exhibit different symmetry breaking compared to that found in previous studies[4,19]. The spin texture in the SkX with does not have the sixfold rotational symmetry in addition to the inversion symmetry, as shown in Fig. 1c, which is different from that in chiral and frustrated magnets in Fig. 1a. We here call this state the threefold-rotational-symmetric SkX (T-SkX). Meanwhile, the SkX shows the inversion symmetry breaking. Furthermore, we elucidate that topological phase transitions between the T-SkX, the SkX, and another non-topological triple-Q state are caused by a change with respect to the relative phase in Eq. (1), which is controlled by the degree of the mirror symmetry breaking. This mechanism for the SkXs might be useful to understand a microscopic origin of the SkX in GdPdSi[36-38], as the underlying lattice structure without the mirror plane on a magnetic layer is common[39].

Results

Model

Let us start by showing an effective spin model starting from a multi-band periodic Anderson model consisting of localized and itinerant electrons. To trace out the itinerant electron degree of freedom and obtain effective magnetic interactions between localized spins, we adopt the standard Schrieffer-Wolff transformation[40] and perform the perturbative expansion of the grand potential with respect to the anisotropic spin-charge coupling[24,41], as detailed in Supplementary Information. Generally, the effective spin model is given bywhere , is the Fourier transform of the localized electron spin at site i (), and the coefficient 2 arises from the contribution. The effective spin model consists of an isotropic RKKY interaction, symmetric anisotropic RKKY interaction, and antisymmetric DM-type RKKY interaction with coupling constants , , and , respectively. The coupling constants are defined bywhere corresponds to the bare susceptibility of itinerant electrons, is the Kronecker delta, and is the Levi-Civita symbol. In Eq. (2), the wave vector is chosen by supposing that , which is relevant to the lattice symmetry. The anisotropic interactions, , and , originate from the atomic spin-orbit coupling[42-44]. The number of and nonzero components of the interactions are determined by the lattice symmetry.

For the above effective spin model, we consider the lattice structure in Fig. 1d consisting of a magnetic layer sandwiched by two nonmagnetic layers. The nonmagnetic ions at () are located above (below) the downward (upward) triangles on the magnetic layer at , which breaks the horizontal mirror symmetry at while keeping the inversion symmetry. The lattice symmetry is compatible with the point group symmetry. In this situation, we set three and four independent coupling constants to satisfy the symmetry. The former is given by , , and and the latter is given by , , , and : , , , , , (all other coupling constants are zero). Among three anisotropic coupling constants, we focus on the effect of , which originates from the horizontal mirror symmetry breaking and is characteristic of the symmetry, on the stabilization of the multiple-Q states, and and are neglected for simplicity[31,34]. In the end, the effective spin model is summarized as

Here, , , ==1, ==, ==, ==, and == (all other component of are zero). The symmetric aniostropic interaction with is qualitatively different from the antisymmetric DM interaction: the former can appear irrespective of the inversion symmetry, while the latter requires the inversion symmetry breaking, and thus vanishes in Eq. (6). The term also appears in the other trigonal crystal systems. We also introduce the Zeeman coupling to an external magnetic field H along the z direction.

Magnetic phase diagram

A magnetic phase diagram of the model in Eq. (6) is calculated by simulated annealing combined with the standard Metropolis local updates. Figure 2a shows the magnetic phase diagram while changing and H in the unit of J at a temperature of 0.01. To identify magnetic phases, we compute the magnetization and the spin structure factor , where is the position vector at site j, is the system size, and is the thermal average. We also calculate the spin scalar chirality where the subscript R represents the center of the triangle and i, j and k are in the counterclockwise order. We obtain six different magnetic phases besides the single-Q (1Q) conical state for and the fully-polarized (FP) state for , whose real-space spin configurations and the in-(out-of-)plane spin structure factor [] are shown in Fig. 2b–g. We also show the skyrmion density for each spin configuration in Fig. 2h–m.

Figure 2

Magnetic phase diagram and characteristics of magnetic phases. (a) -H magnetic phase diagram for the model in Eq. (6) in the unit of J. The 3Q-I, 3Q-II, 3Q-III, SkX-2, SkX-1, 3Q-Ch, 1Q conical, and FP represent the triple-Q I, triple-Q II, triple-Q III, SkX, T-SkX, triple-Q chiral, single-Q conical, and fully polarized states, respectively. In the hatched region, energies for several magnetic states are degenerate and it is difficult to determine the phase boundaries. (b–g) Snapshots of the spin configurations in (b) 3Q-I for and , (c) 3Q-II for and , (d) SkX-2 for and , (e) SkX-1 for and , (f) 3Q-III for and , and (g) 3Q-Ch for and . The arrows and contour denote the xy and z components of the spin moments, respectively. The square root of in-plane and out-of-plane spin structure factors in the Brillouin zone are shown in upper and lower panels, respectively, where the dashed circles highlight , , and and the component is removed for better visibility. (h–m) Real-space distributions of the skyrmion density for the spin configurations in (b–g), respectively. (n) H dependences of the magnetization (red square) and spin scalar chirality (blue circle) for (filled symbols) and (open symbols).

For , the model in Eq. (6) reduces to the isotropic RKKY model, which stabilizes the 1Q conical state for any H. By introducing , the multiple-Q instabilities occur: the triple-Q I (3Q-I) state is stabilized for small H, while the triple-Q II (3Q-II) state is stabilized for large H, as shown in Fig. 2a. Their spin modulations are mainly characterized by the in-plane single-q component, which smoothly connects to the 1Q conical state. Meanwhile, they exhibit different peak structures in , as shown in Fig. 2b,c: there is a dominant peak at in the 3Q-I state, whereas there are two dominant peaks at and in the 3Q-II state in addition to the peak at in . Both phases are topologically trivial without .

While increasing , the 3Q-I state is replaced by the SkX (SkX-2) in the low-field region for and the T-SkX (SkX-1) in the intermediate-field region for , as shown in Fig. 2a. Both SkXs are characterized by the triple-Q peaks with the same intensities, as shown in Fig. 2d,e. By looking into the real-space spin configurations, they are formed by a vortex with vorticity and two vortices with in a magnetic unit cell in Fig. 2d,e, which indicates the inversion symmetry breaking. The positions at the cores with negative are different with each other. They are located at the cores with () for the SkX ( T-SkX). Such a difference results in the different skyrmion numbers, which is clearly found in Fig. 2j,k.

In the high-field region, the 3Q-II state and T-SkX are replaced by the other topologically trivial triple-Q states depending on : the triple-Q III (3Q-III) or the triple-Q chiral (3Q-Ch) state. The 3Q-III state is mainly characterized by the in-plane double-Q peaks in Fig. 2f, while the 3Q-Ch state is by the in-plane triple-Q peaks with equal intensities in Fig. 2g. The 3Q-Ch state exhibits nonzero , although the skyrmion number becomes zero.

We show H dependences of M and for and 0.2 in Fig. 2n. While increasing H, jumps of M and appear when the skyrmion number changes: Two jumps between the T-SkX and the 3Q-I are found at and for and other two jumps are found between the SkX and the T-SkX and between the T-SkX and the 3Q-Ch for . The transition between the T-SkX and the 3Q-III for also shows jumps of M and .

Mechanism of the topological transition

Next, we show the transformation of the skyrmion numbers on the basis of the phase degrees of freedom among the constituent triple-Q density waves. We find that the spin configurations for the T-SkX in Fig. 2e and SkX in Fig. 2d are summarized in a single expression aswhere A and are additional variational parameters compared to Eq. (1), and we set and for . There are two types of phase degrees of freedom in Eq. (7). One is a phase among the constituent waves, , which induces the transformation between the SkX and the vortex crystal with staggered spin scalar chirality[45]. The other is a relative phase between in-plane- and z-spin components, , which induces the different types of the SkX, as discussed in the introduction: the T-SkX for and SkX for , where depends on the other variational parameters. From the symmetry viewpoint, nonzero in the SkX and T-SkX shows the sixfold-rotational symmetry breaking in addition to the inversion symmetry breaking, which is in contrast to in the SkX in chiral magnets.

The phase diagram obtained by the simulated annealing in Fig. 2a is well reproduced by the variational spin ansatz in Eq. (7) for . This means that is evaluated through the spin ansatz in Eq. (7). Figure 3a shows H dependence of by simulated annealing and variational calculation at . Both the results show that the SkX exhibits , the T-SkX exhibits , and the 3Q-Ch exhibits . Thus, the phase transitions among the SkX, the T-SkX, and the 3Q-Ch are regarded as the phase transitions with respect to . In other words, the topological transitions are caused by the changes of the types of constitute waves, as shown in Fig. 3a.

Figure 3

Mechanism of the topological transition for . (a) H dependence of the phase obtained from the simulated annealing (SA) and variational calculation (VC) in the model in Eq. (6) at . Schematic pictures of the constitute waves in each SkX are shown in the inset. (b) dependences of energies from different contributions, , , and in the unit of J, , and H, respectively, for , , , , and in Eq. (7) with . The vertical dashed line represents the boundary between the SkXs with and 2.

The phase transition characterized by the change of is due to the mirror symmetry breaking in the lattice structure. To demonstrate that, we calculate the energy contributions from each term in the model in Eq. (6): the RKKY energy , the anisotropic energy , and the Zeeman energy . Figure 3b shows dependences of , , and for the spin ansatz in Eq. (7) with at fixed , , , , and , where is around . has little dependence, whereas and show distinct behaviors against ; () decreases while increasing (decreasing) . In other words, arising from the mirror symmetry breaking tends to favor the SkX with , while the magnetic field tends to favor the SkX with . The competition between these distinct behaviors causes the filed-induced transition from the SkX to the T-SkX with .

Summary

In conclusion, we clarify that the magnetic anisotropy arising from the breaking of the mirror symmetry is another way to stabilize the SkXs in itinerant magnets irrespective of the spatial inversion symmetry. On the basis of simulated annealing and variational calculation, we show that the anisotropic RKKY interaction induces two SkXs with different topological numbers, which accompanies the spontaneous inversion symmetry breaking. Moreover, we find that two SkXs are transformed with each other by changing the anisotropic RKKY interaction and magnetic field, the former of which is tuned by the degree of mirror symmetry breaking.

Our study reveals that the and SkXs are stabilized even without the multi-spin interaction in itinerant magnets, which is distinct from the previous one in the Kondo lattice model without the magnetic anisotropy[23]: The anisotropic bilinear exchange interaction plays an important role in the stabilization of the former SkXs, while the isotropic biquadratic interaction is important for the latter one[24]. Although both the systems exhibit similar skyrmion textures, the degeneracy in terms of the vorticity and helicity is different owing to the different mechanisms. The SkXs by the isotropic biquadratic interaction are energetically degenerate for different vorticity and helicity, while the present SkXs have a definite vorticity and helicity depending on the sign of the anisotropic interaction. Reflecting such a difference, the SkXs in the present model induce a different Goldstone mode from that in the previous model, which results in different dynamics. Furthermore, the anisotropic response against the electromagnetic field is anticipated due to the nature of magnetic anisotropy, which might give rise to further unconventional multiple-Q states. Such a theoretical exploration of the SkXs based on magnetic anisotropy will be left for future study.

Finally, let us discuss a relevant material in the present mechanism. The centrosymmetric itinerant magnet GdPdSi[36] might be a candidate material, which hosts the skyrmion crystal in an external magnetic field. The importance of the RKKY interaction from the nesting of the Fermi surfaces has already been suggested by the angle-resolved photoemission spectroscopy[46,47]. In addition, the anisotropic RKKY interaction would play an important role in this compound, as the magnetic Gd ions form the triangular lattice and they are sandwiched by the nonmagnetic Pd and Si so that the mirror symmetry on each magnetic layer is broken[39,48]. Indeed, the importance of the multi-orbital degrees of freedom, which become the microscopic origin of the anisotropic RKKY interaction, has also been implied by first-principle calculations[46,49]. It would be interesting to test our scenario for the SkX in GdPdSi by considering the superstructure of the Pd and Si and the effect of the spin-orbit coupling. Our mechanism will also shed light on engineering the SkXs in quasi-two-dimensional magnetic materials including surface, domain, and layered systems.

Methods

Simulated annealing

We perform the simulated annealing combined with the standard Metropolis local updates under the periodic boundary condition. In the simulation, we gradually reduce the temperature with a rate , where is the temperature at the nth step. We set the initial temperature and the coefficient . The final temperature is reached after total steps, where we perform Monte Carlo sweeps at each temperature. At the final temperature, we perform Monte Carlo sweeps for thermalization and measurements, respectively. Although Figs. 2 and 3a show the results for , we confirm that the obtained result does not change for . We also confirm that the simulations with different values of , , and 0.95499, give the same result.

Variational calculation

We here present the details of the spin ansatz in Eq. (7) and the variational calculation. To find the spin configuration in Eq. (7), we start from a general spin ansatz of the single-Q spiral state given bywhere , , and (). This spin configuration expresses an elliptical wave, where the axis with a length of 2a is parallel to and the axis with a length of 2b () is parallel to . The spin ansatz in Eq. (8) describes a variety of spin density waves depending on and : the spiral wave (), standard elliptical wave ( and ), rotated elliptical wave ( and ), standard sinusoidal wave ( and ), and rotated sinusoidal wave ( and ), where “standard” means that the axes are parallel to or . A schematic picture of spiral plane in Eq. (8) is shown in Supplementary Information.

We rewrite the spin configuration in Eq. (8) aswhere , , and satisfy , , , and , respectively. Then, the spin ansatz in Eq. (7) is obtained by superposing three spin density waves in Eq. (9) in addition to a uniform z component, which is given by

The variational parameters A, , , and in Eq. (7) are related to that in Eq. (10) as , and . From the definition, does not depend on , i.e., .

In the variational calculations in Fig. 3a, we optimize a, b, , , and as the variational parameters for . After obtaining the optimal parameters, we calculate in Fig. 3a from the difference between phases of and . In Supplementary Information, we show H dependences of the magnetic moment, magnetization, and spin scalar chirality by the variational calculation and the simulated annealing. Compared to the results, one can find that the variational spin ansatz in Eq. (7) corresponds to the spin textures obtained by the simulated annealing.

Supplementary Information.