Temperature-driven transition from skyrmion to bubble crystals in centrosymmetric itinerant magnets

Interplay between itinerant electrons and localized spins in itinerant magnets gives rise to a variety of noncoplanar multiple-$Q$ spin textures, such as the skyrmion, hedgehog, meron, and vortex. We elucidate that another type of multiple-$Q$ state consisting of collinear sinusoidal waves, a magnetic bubble crystal, appears at finite temperatures in a centrosymmetric itinerant electron system. The results are obtained for the classical Kondo lattice model with easy-axis single-ion anisotropy on a triangular lattice by a large-scale numerical simulation. We find that a finite-temperature topological phase transition between the skyrmion crystal and the bubble crystal occurs by changing the temperature. We obtain the minimal key ingredients for inducing the finite-temperature transition by analyzing an effective spin model where it is shown that the synergy between the multiple-spin interaction and magnetic anisotropy plays a significant role.


Introduction
A magnetic skyrmion, which is characterized by a swirling spin texture to have a nontrivial topological (skyrmion) number, has been extensively studied in condensed matter physics, as it exhibits a number of intriguing physical phenomena that arises from emergent magnetic fields, such as the topological Hall effect [1,2,3,4]. Since the discovery of the periodic array of the magnetic skyrmion dubbed a skyrmion crystal (SkX) in MnSi in 2009 [5,6], various types of the SkXs have been observed in metals [5,7,8,9,10], semiconductors [8,11], and insulators [12,13,14,15,16] in noncentrosymmetric magnets [17]. Furthermore, a fertile playground for the SkXs has been established in recent years by detecting them in centrosymmetric magnets [18,19,20,21,22]. As the topological properties in the swirling spin texture are robust against external stimuli, the SkXs are promising for high-efficient storage and memory applications [23,24,25,26,27,28].
Theoretically, the stabilization mechanisms of the SkXs have also been clarified in various electron and/or spin systems. The most intuitive model to stabilize the SkXs is the localized spin model consisting of the ferromagnetic exchange interaction and the Dzyaloshinskii-Moriya (DM) interaction, the latter of which arises from the relativistic spin-orbit coupling without the inversion symmetry [29,30]. The model at a zero field exhibits a magnetic spiral state originating from the competition between the ferromagnetic and DM interactions where the ferromagnetic (DM) interaction tends to align (tilt) the neighboring spins. The system undergoes a magnetic field-induced phase transition to the SkXs, which are formed by superposing three spiral states [31,32]. Subsequently, it is shown that the SkXs are also stabilized by the other interactions and fluctuations even without the Dzyaloshinskii-Moriya interaction [33]. For example, the spin model with the frustrated exchange interactions and/or magnetic anisotropy gives rise to the SkXs [34,35,36,37,38]. Another example is found in the itinerant electron model which implicitly includes effective multiple-spin interactions arising from the Fermi surface instability [39,40,41,42,43]. Notably, these mechanisms enable the formation of the SkXs with short-period magnetic modulations in contrast to that induced by the DM mechanism where the magnetic modulation usually has a long period owing to the small DM interaction compared to the exchange interaction. The smallperiod SkX might be important for next-generation spintronic applications, as it realizes an energy-efficient information storage based on high density skyrmions [28].
In the present study, we investigate the stability of the SkX in centrosymmetric itinerant magnets at finite temperatures. By considering the classical Kondo lattice model with an easy-axis single-ion anisotropy on a triangular lattice and performing large-scale Langevin dynamics simulations combined with the kernel polynomial method (KPM-LD) [101], we find that the system exhibits a finite-temperature phase transition between the SkX described by the triple-Q spiral waves and the bubble crystal described by the triple-Q collinear sinusoidal waves. This phase transition is regarded as the topological phase transition where the scalar chirality is turned on and off. We show that such a finite-temperature phase transition between the SkX and the bubble crystal is captured by an effective spin model derived from the Kondo lattice model, suggseting that the interplay between the multiple-spin interactions and magnetic anisotropy is essential to the transition. Our results indicate that temperature fluctuations can be a source of multiple-Q states and drive a further topological phase transition in the skyrmion-hosting itinerant magnets.

Results and discussion
We discuss the stability of the SkX in itinerant magnets by performing large-scale numerical simulations. We show the result in the Kondo lattice model on a triangular lattice in section 2.1 and in the effective spin model in section 2.2.

Classical Kondo lattice model
To analyze the finite-temperature effect on the SkX in centrosymmetric itinerant magnets, we consider the classical Kondo lattice model (or the double exchange model) on the triangular lattice, where the previous studies have shown that the SkX is stabilized at zero temperature [43,102,103]. The Kondo lattice model consists of itinerant electrons and localized spins, whose Hamiltonian is given by where c † iσ (c iσ ) is a creation (annihilation) operator of an itinerant electron at site i and spin σ. The first term in (1) stands for the kinetic energy of itinerant electrons where t ij is a hopping parameter between sites i and j. The second term in (1) stands for the exchange coupling between itinerant electron spins where J is a coupling constant and σ is the vector of Pauli matrices. We consider the localized spin S i as the classical one with |S i | = 1 for simplicity. The third term in (1) represents the easy-axis single-ion anisotropy A > 0. The fourth term in (1) represents the Zeeman coupling to localized spins under an external magnetic field along the z direction.
The model in (1) shows a plethora of multiple-Q states including the SkXs by changing t ij , J, A, H, and the chemical potential µ [39,104,43,102]. Among them, we focus on a parameter set where two SkXs are robustly stabilized in the ground state by choosing the nearest-neighbor hopping t 1 = 1, the third nearest-neighbor hopping t 3 = −0.85, J = 1, and µ = −3.5; the SkX with the skyrmion number of two (one) appears at the zero (finite) magnetic field for A = 0 [43,103]. The emergent SkXs are characterized by the superposition of the triple-Q spin density waves with Q 1 = (π/3, 0), Q 2 = (−π/6, √ 3π/6), and Q 3 = (−π/6, − √ 3π/6), which are dictated by the peaks of the bare susceptibility of itinerant electrons. By introducing A, the stability region of the SkX with the skyrmion number of one under the magnetic fiend is considerably extended compared to that with the number of two [102]. We here examine the temperature effect on the SkX with the skyrmion number of one in the presence of A, since all the observed SkXs in itinerant magnets correspond to the SkX with the number of one [19,20,21]. The temperature effect on the SkXs for A = 0 has been investigated by the author and his collaborators where it was shown that the SkX changes into the topologically-trivial triple-Q state breaking the threefold rotational symmetry at finite temperatures [55].
The behavior of the finite temperature is investigated by performing the KPM-LD simulation, which is one of the efficient methods in the itinerant electron model [105,101]. This method is used for the large-scale system where the partition function consists of itinerant electrons coupled to classical boson fields, which has been applied to similar itinerant electron models [101,106,71,43,107,108,109,110]. In the present Kondo lattice Hamiltonian with the classical localized spins in (1) is the free energy for a given localized spin configuration and is given by Here, ρ(ω) is the density of states of H under {S i }, which is evaluated by the kernel polynomial method [105]. We expand the density of states by up to 2000th order of Chebyshev polynomials with 16 2 random vectors [111]. Meanwhile, Tr {S i } is evaluated by the spin configuration generated by the Langevin equation, where we adopt a projected Heun scheme [112] for 1000-5000 steps with the time interval ∆τ = 2. We start the simulation from initial states with random spin configurations for the system sizes with N = L 2 at L = 96 and 120 with periodic boundary conditions in both directions. The results at different temperatures are obtained independently starting from different random spin configurations. Figure 1 shows the result at A = 0.007 and H = 0.008 obtained by the KPM-LD simulations for the system size with L = 96 and 120. To identify the magnetic phases, we compute the magnetic moments with wave vector q, where S α (q) is the α component of the spin structure factor. We also calculate the spin scalar chirality as where µ = (u, d) represent upward and downward triangles, respectively, and R and R ′ represent the position vectors at the centers of triangles; χ R = S j · (S k × S l ) is the local spin chirality at R, where j, k, l are the sites located at the triangle vertices at R in the anticlockwise order. The nonzero uniform component of χ q , i.e., χ 0 , can indicate the magnetic phase with the topologically nontrivial spin texture, such as the SkX. It is noted that the long-range order in terms of the xy spin is limited to zero temperature in the present two-dimensional system owing to the Mermin-Wagner theorem [113], while the z spin and scalar chirality can order at finite temperatures. As shown in figure 1, the system has nonzero χ 0 at zero temperature. While increasing temperature, χ 0 gradually decreases and vanishes around T 1 ∼ 0.0045, which indicates the topological phase transition regarding the spin scalar chirality degree of freedom. In other words, this phase transition means the spontaneous Z 2 mirror symmetry breaking of the Kondo lattice Hamiltonian according to the ordering of the scalar chirality. Meanwhile, the spin momentm z Q takes a finite value above T 1 wherẽ m z Q is defined as where Q 1 -Q 3 represents the ordering vectors that the bare susceptibility of itinerant electrons is maximized. Although the magnetic state is characterized by the state with the multiple peaks at Q ν in the case of A = 0 [43], the peak positions can be slightly moved from Q ν when considering the effect of A and T . Indeed, such situations have been found in the frustrated magnets with the strong easy-axis anisotropy [37] and the axial next-nearest-neighbor Ising model where the anisotropy is infinitesimally large [114,115,116]. To take into account such a change of the peak positions, we introduce q n νη in (6) for the nth-neighbor displacement vectors measured from Q ν (η is the number of the nth-neighbor displacement vectors), e.g., η = 0 for n = 0 and η = 1-6 for n = 1-3. We here take the summation up to n = 3. The transition temperature wherem z Q vanishes is around T 2 ∼ 0.009, which is clearly higher than T 1 . This indicates that the phase transition at T 2 means the ordering of the z spin component caused by the single-ion anisotropy. Thus, the transitions in terms of the spin and chirality occur at different temperatures in the model in (1). Such a separation between spin and chirality degrees of freedom has been discussed in various itinerant electron models, such as the Hubbard model on the kagome lattice [117,118] and the Kondo lattice model without the magnetic anisotropy on the triangular lattice [119,55]. The reason why T 1 and T 2 are different might be attributed to the presence of two different energy scales arising from the effective magnetic interaction mediated by itinerant electrons and the single-ion anisotropy. Indeed, similar phase transitions are found in the effective spin model derived to include these two factors, as discussed in section 2.2.
To examine the spin textures in each phase, we plot the real-space spin and chirality configurations obtained by the KPM-LD simulations in figure 2. To integrate out the short-wavelength fluctuations, we average over 200-400 spin configurations for different time steps. In the low-temperature phase at T = 0.002, the SkX spin texture is found in figure 2(a); the skyrmion cores at S z i ∼ −1 form a triangular lattice with the lattice constant 4π/( √ 3|Q ν |), which indicates the triple-Q peak structures at q ≃ Q 1 , Q 2 , and Q 3 with equal intensities. This phase also accompanies a uniform scalar chirality, as shown in figure 2(b), where the skyrmion number becomes ±1. This SkX has the degeneracy with respect to the helicity and vorticity, which means that this state can take both positive and negative chirality depending on the initial spin configuration, which is similar to that in frustrated magnets [34,35,36,37]. It is noted that the degeneracy between the SkXs with positive and negative scalar chirality is lifted if the system has the Dzyaloshinskii-Moriya interaction [120] and bond-dependent anisotropic interaction [38,121,92,62,56,122] that might not be neglected in real materials.
Figures 2(c) and 2(d) represent the spin and chirality distributions in the intermediate temperature region, respectively, wherem z Q is nonzero but χ 0 vanishes. The spin configuration in figure 2(c) shows the almost collinear-type spin texture, which results in no scalar chirality in figure 2(d). Nevertheless, the z component of the spin moment still forms the triangular lattice, which is clearly found in the inset of figure 2(d). Thus, the triple-Q collinear spin texture is realized in the intermediate region, which corresponds to the magnetic bubble crystal. The bubble core disappears in the paramagnetic state above T 2 , as shown in figures 2(e) and 2(f).
The above result indicates that the finite-temperature phase transition between the SkX and the bubble crystal is characterized by nonzero χ 0 . Both the states are characterized by the triple-Q spin density waves, although the type of the individual spin density waves is different with each other. On the one hand, the SkX consists of the three spirals whose plane is lied along the xz or yz plane, which is given by where e 1 =x, e 2 = −x/2 + √ 3ŷ/2, and e 3 = −x/2 − √ 3ŷ/2 (x andŷ are the unit vectors along the x and y directions). I xy , I z , m z , and θ ν depends on the model parameters. On the other hand, the bubble crystal is described by the three sinusoidal waves which oscillates along the z direction as Thus, one can find that the SkX spin texture reduces to the bubble one in the limit of I xy → 0, which indicates that the scalar chirality is largely suppressed in the bubble crystal.
Finally, let us comment on numerical simulations. Although we have adopted the efficient KPM-LD simulations, which roughly reduce the computational cost of the conventional Monte Carlo sampling with the exact diagonalization by N 3 , it still needs a huge computational cost for the large system sizes with N ∼ 10 4 . Depending on the initial spin configurations, the spin states are easily trapped by the multi-domain structures. In order to show the variance for the different initial spin configurations at the same temperature, we show the averaged values ofm z Q and χ 0 denoted asm z,ave Q and χ ave 0 , respectively, by performing seventeen independent simulations in the SkX phase at T = 0.002 and in the bubble crystal phase at T = 0.006, as shown in figure 1. The result at T = 0.006 shows a good convergent behavior, while that at T = 0.002 shows a large variance, especially for χ ave 0 . This is attributed to the spin and/or chiral domain structures, the latter of which arises from the degeneracy between the SkX with the negative and positive chiralities. Thus, although a clear signal of the transition between the SkX and bubble crystal is found in the present study, a more careful finite-size scaling is required to examine the nature of the phase transition, which will be left for future study.

Effective spin model
To investigate key ingredients for the finite-temperature transition between the SkX and bubble crystal, we investigate an effective spin model in the Kondo lattice model. The effective spin model is derived from the perturbation expansion in terms of J [42], which is given by Temperature-driven transition from skyrmion to bubble crystals in centrosymmetric itinerant magnets9 where the term in the square bracket corresponds to the effective interaction arising from the spin-charge coupling in the Kondo lattice model, while the second and third terms are the same as the third and fourth terms in (1). The first term represents the spin interactions in momentum space;J is the second-order contribution with respect to J, while K is the fourth-order contribution with respect to J, both of which are positive. We consider the specific q components of the spin interactions in momentum space where S q = (1/ √ N ) i S i e −iq·r i is the Fourier component of S i ; we only take into account the interactions at Q 1 , Q 2 , and Q 3 , which are dictated by the peaks of the bare susceptibility of itinerant electrons in the original Kondo lattice model. This model is contrast to a spin model with similar biquadratic interactions in real space, which also describes the multiple-Q states [123,9,124,120]. We setJ = 1 as the energy unit of the model in (9). The other parameters are chosen at K = 0.5, A = 0.3, and H = 0.7 so as to stabilize the SkX with the skyrmion number of one at zero temperature, and the system size is taken at L = 96.
The following result is calculated by performing simulated annealing from high temperature, which is based on the standard Metropolis local updates in real space. Similar to the previous studies [42,125,59], we gradually reduce the temperature with the rate α = 0.99995-0.99999 step by step, starting from the initial temperature T 0 = 1.5 for a random spin configuration. At the target temperature, we perform 10 5 -10 6 Monte Carlo sweeps for measurements after equilibration. The following data at different temperature are obtained for the different initial spin configurations. The calculation of the effective spin model can be performed by a considerably smaller computational cost than that by the KPM-LD simulation of the Kondo lattice model.  the single domain in each magnetic phase, which might be due to the interactions limited at particular Q ν . The phase sequence while changing the temperature is similar to that in the Kondo lattice model in figure 1. The SkX is stabilized at low temperatures, the bubble crystal is stabilized at intermediate temperatures, and the paramagnetic state appears at high temperatures. We show the real-space snapshots of the spin and chirality in each phase in figures 4(a)-4(f), which well correspond to those in the Kondo lattice model. In other words, the separation of the transition temperature in terms of spin and chirality degrees of freedom also arises in the effective spin model, which indicates that the effective spin model qualitatively captures the finite-temperature behavior of the Kondo lattice model. Thus, the model in (9) is a minimal model to describe the finite-temperature transition from the SkX to the bubble crystal. In fact, we did not obtain such a phase transition when we omit any ofJ, K, A, and H.

Summary
To summarize, we have investigated the stability of the SkX in the Kondo lattice model with the easy-axis single-ion anisotropy on the centrosymmetric triangular lattice. By performing KPM-LD simulations for the large system size, we find that thermal fluctuations drive the SkX to the bubble crystal. Although both states are characterized by the triple-Q states where the cores form the triangular lattice, only the SkX has the net scalar chirality inducing the topological Hall effect, while the bubble crystal has a collinear spin texture. We also show that finite-temperature transition between the SkX and the bubble crystal is found in the effective spin model derived from the Kondo lattice model, which suggests that the interplay between the biquadratic interaction and the single-ion anisotropy plays an important role in the transition. As the similar phase transition is also found in frustrated magnets [37], the transition between the SkX and the bubble crystal is a universal feature in centrosymmetric magnets with the strong magnetic anisotropy.