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\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_2$$\end{document}2PdSi\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$_3$$\end{document}3, but also for exploring further skyrmion-hosting materials in trigonal itinerant magnets.

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][5][6] . The SkX exhibits a nonzero topological winding number called the skyrmion number N sk , which is defined as N sk = R � R /4π , where R is a skyrmion density related to the solid angle consisting of three spins S i , S j , and S k on the triangle R: tan (� R /2) = S i · (S j × S k )/(1 + S i · S j + S j · S k + S k · S i ) 7 . The study of the SkX has attracted much attention, as the swirling topological magnetic texture owing to nonzero N sk gives rise to an emergent electromagnetic field through the spin Berry phase and results in intriguing transport phenomena and dynamics [8][9][10][11][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) as where e η and e z are the unit vectors along the in-plane and z directions, respectively. Q ηi = Q η · r i + φ η , and Q ′ ηi = Q ηi + ψ η 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 e η � e z × Q η ( e η Q η ) and ψ 1 = ψ 2 = ψ 3 = 0 or π represents the Bloch-type (Néel-type) SkX, while that for ψ 1 = ψ 2 = 0 and ψ 3 = π represents the anti-type SkX. The realspace spin texture for the Bloch-type SkX is shown in Fig. 1a. All the SkXs have the skyrmion number of one, n sk ≡ |N sk | = 1 , in the magnetic unit cell and breaks the spatial inversion symmetry irrespective of e η and ψ η 8 . We call them the n sk = 1 SkXs. The n sk = 1 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][21][22] .
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 n sk = 1 SkX, this spin texture exhibits the skyrmion number of two in a magnetic unit cell ( n sk = 2 SkX), whose spatial inversion and/or sixfold rotational symmetries are broken depending on φ η on a discrete lattice. For example, the n sk = 2 SkX with φ η = π shown in Fig. 1b has the inversion symmetry, but the n sk = 2 SkX with φ 1 = 4π/3 , φ 2 = 2π/3 , and φ 3 = π shows the inversion symmetry breaking. Although the n sk = 2 SkX seems to be rare compared to the n sk = 1 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 n sk = 2 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][29][30] arising from the absence of the mirror symmetry on a magnetic layer [ Fig. 1d] induces the SkXs with n sk = 1 and n sk = 2 . The anisotropic RKKY interaction stabilizes the SkXs even without the DM, competing exchange, and multi-spin interactions 25,[31][32][33][34][35] . The obtained SkXs exhibit different symmetry breaking compared to that found in previous studies 4,19 . The spin texture in the SkX with n sk = 1 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 n sk = 1 threefold-rotational-symmetric SkX (T-SkX). Meanwhile, the n sk = 2 SkX shows the inversion symmetry breaking. Furthermore, we elucidate that topological phase transitions between the n sk = 1 T-SkX, the n sk = 2 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 Gd 2 PdSi 3 36-38 , as the underlying lattice structure without the mirror plane on a magnetic layer is common 39 .

Results
Model. Let   www.nature.com/scientificreports/ 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 by where α, β = x, y, z , S Q η is the Fourier transform of the localized electron spin S i at site i ( |S i | = 1 ), and the coefficient 2 arises from the −Q η contribution. The effective spin model consists of an isotropic RKKY interaction, symmetric anisotropic RKKY interaction, and antisymmetric DM-type RKKY interaction with coupling constants J Q η , K αβ Q η , and D α Q η , respectively. The coupling constants are defined by where χ αβ q corresponds to the bare susceptibility of itinerant electrons, δ αβ is the Kronecker delta, and ǫ αβγ is the Levi-Civita symbol. In Eq. (2), the wave vector Q η is chosen by supposing that χ Q η > χ q , which is relevant to the lattice symmetry. The anisotropic interactions, K αβ Q η , and D α Q η , originate from the atomic spin-orbit coupling [42][43][44] . The number of Q η 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 z = c ( z = −c ) are located above (below) the downward (upward) triangles on the magnetic layer at z = 0 , which breaks the horizontal mirror symmetry at z = 0 while keeping the inversion symmetry. The lattice symmetry is compatible with the D 3d point group symmetry. In this situation, we set three Q η and four independent coupling constants to satisfy the D 3d symmetry. The former is given by Q 1 = (2π/6, 0, 0) , Q 2 = (−π/6, √ 3π/6, 0) , and Q 3 = (−π/6, − √ 3π/6, 0) and the latter is given by (all other coupling constants are zero). Among three anisotropic coupling constants, we focus on the effect of K yz Q 1 , which originates from the horizontal mirror symmetry breaking and is characteristic of the D 3d symmetry, on the stabilization of the multiple-Q states, and K xx Q 1 and K yy Q 1 are neglected for simplicity 31,34 . In the end, the effective spin model is summarized as 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 M = (1/N) j �S z j � and the spin structure factor S αα where r j is the position vector at site j, N = 48 2 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 Ŵ = 0 and the fully-polarized (FP) state for H 2 , whose real-space spin configurations and the in-(out-of-)plane spin structure factor S ⊥ s (q) = S xx s (q) + S yy s (q) [ S zz s (q) ] are shown in Fig. 2b-g. We also show the skyrmion density R for each spin configuration in Fig. 2h-m.
For Ŵ = 0 , 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. www.nature.com/scientificreports/ Meanwhile, they exhibit different peak structures in S zz s (q) , as shown in Fig. 2b,c: there is a dominant peak at Q 2 in the 3Q-I state, whereas there are two dominant peaks at Q 1 and Q 3 in the 3Q-II state in addition to the peak at Q 2 in S ⊥ s (q) . Both phases are topologically trivial without χ sc . While increasing Ŵ , the 3Q-I state is replaced by the n sk = 2 SkX (SkX-2) in the low-field region for Ŵ 0.1 and the n sk = 1 T-SkX (SkX-1) in the intermediate-field region for Ŵ 0.05 , 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 realspace spin configurations, they are formed by a vortex with vorticity ν = −2 and two vortices with ν = +1 in a magnetic unit cell in Fig. 2d,e, which indicates the inversion symmetry breaking. The positions at the cores with negative S z i are different with each other. They are located at the cores with ν = −2 ( ν = +1 ) for the n sk = 2 SkX ( n sk = 1 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 n sk = 1 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 χ sc , although the skyrmion number becomes zero.
We show H dependences of M and |χ sc | for Ŵ = 0.075 and 0.2 in Fig. 2n. While increasing H, jumps of M and χ sc appear when the skyrmion number changes: Two jumps between the n sk = 1 T-SkX and the 3Q-I are found

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 n sk = 1 T-SkX in Fig. 2e and n sk = 2 SkX in Fig. 2d are summarized in a single expression as where A and m z are additional variational parameters compared to Eq. (1), and we set e η � e z × Q η and ψ ≡ ψ 1 = ψ 2 = ψ 3 for Ŵ > 0 . 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 n sk = 1 T-SkX for 0 < ψ < ψ c and n sk = 2 SkX for ψ c ≤ ψ ≤ π/2 , where ψ c depends on the other variational parameters. From the symmetry viewpoint, nonzero ψ in the n sk = 2 SkX and n sk = 1 T-SkX shows the sixfold-rotational symmetry breaking in addition to the inversion symmetry breaking, which is in contrast to ψ = 0 in the n sk = 1 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 Ŵ 0.2 . 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 Ŵ = 0.2 . Both the results show that the n sk = 2 SkX exhibits ψ = 0.5π , the n sk = 1 T-SkX exhibits 0.24π ψ 0.3π , and the 3Q-Ch exhibits ψ = 0.5π . Thus, the phase transitions among the n sk = 2 SkX, the n sk = 1 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.
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 E RKKY = −2J η S Q η · S −Q η /N , the anisotropic energy E Ŵ = −2Ŵ ηαβ I αβ Q η S α Q η S β −Q η /N , and the Zeeman energy E Zeeman = −H i S z i /N . Figure 3b shows ψ dependences of E RKKY , E Ŵ , and E Zeeman for the spin ansatz in Eq. (7) with |S i | = 1 at fixed A = 1/ √ 2 , m z = 0 , φ 1 = 4π/3 , φ 2 = 2π/3 , and φ 3 = π , where ψ c is around π/3 . E RKKY has little ψ dependence, whereas E Ŵ and E Zeeman show distinct behaviors against ψ ; E Ŵ ( E Zeeman ) decreases while increasing (decreasing) ψ . In other words, Ŵ arising from the mirror symmetry breaking tends to favor the n sk = 2 SkX with ψ = π/2 , while the magnetic field tends to favor the n sk = 1 SkX with ψ = 0 . The competition between these distinct behaviors causes the filed-induced transition from the n sk = 2 SkX to the n sk = 1 T-SkX with 0.24π ψ 0.3π.

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 n sk = 1 and n sk = 2 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 Gd 2 PdSi 3 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 Gd 2 PdSi 3 by considering the superstructure of the Pd and Si and the effect of the spinorbit 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 T n+1 = αT n , where T n is the temperature at the nth step. We set the initial temperature T 0 = 1 and the coefficient α = 0.99954 . The final temperature T ≃ 0.01 is reached after total 10 6 steps, where we perform 10 2 Monte Carlo sweeps at each temperature. At the final temperature, we perform 10 6 Monte Carlo sweeps for thermalization and measurements, respectively. Although Figs. 2 and 3a show the results for N = 48 2 , we confirm that the obtained result does not change for N = 96 2 . We also confirm that the simulations with different values of α , α = 0.99908, 0.99541 , 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 by where Q ′′ ηi = Q η · r i + ϕ η , ẽ ηy = e η cos θ − e z sin θ , and ẽ ηz = e η sin θ + e z cos θ ( 0 ≤ θ < π ). This spin configuration expresses an elliptical wave, where the axis with a length of 2a is parallel to ẽ ηz and the axis with a length of 2b ( a > b ≥ 0 ) is parallel to ẽ ηy . The spin ansatz in Eq. (8) describes a variety of spin density waves depending on r ≡ b/a and θ : the spiral wave ( r = 1 ), standard elliptical wave ( 0 < r < 1 and θ = 0, π/2 ), rotated elliptical wave ( 0 < r < 1 and θ = 0, π/2 ), standard sinusoidal wave ( r = 0 and θ = 0, π/2 ), and rotated sinusoidal wave ( r = 0 and θ = 0, π/2 ), where "standard" means that the axes are parallel to e η or e z . A schematic picture of spiral plane in Eq. (8) is shown in Supplementary Information.
S η i = e η A ⊥ sin(Q ′′ ηi + � ⊥ ) + e z A z cos(Q ′′ ηi + � z ), www.nature.com/scientificreports/ In the variational calculations in Fig. 3a, we optimize a, b, θ , ϕ η , and M z as the variational parameters for N = 12 2 . After obtaining the optimal parameters, we calculate ψ in Fig. 3a from the difference between phases of S y Q 1 and S z Q 1 . 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.

Data availibility
The datasets generated during and/or analysed during the current study are available from the corresponding author on reasonable request.