Locking of skyrmion cores on a centrosymmetric discrete lattice: onsite versus offsite

A magnetic skyrmion crystal (SkX) with a swirling spin configuration, which is one of topological spin crystals as a consequence of an interference between multiple spin density waves, shows a variety of noncoplanar spin patterns depending on a way of superposing the waves. By focusing on a phase degree of freedom among the constituent waves in the SkX, we theoretically investigate a position of the skyrmion core on a discrete lattice, which is relevant with the symmetry of the SkX. The results are obtained for the double exchange (classical Kondo lattice) model on a discrete triangular lattice by the variational calculations. We find that the skyrmion cores in both two SkXs with the skyrmion number of one and two are locked at the interstitial site on the triangular lattice, while it is located at the onsite by introducing a relatively large easy-axis single-ion anisotropy. The variational parameters and the resultant Fermi surfaces in each SkX spin texture are also discussed. The different symmetry of the Fermi surfaces depending on the core position is obtained when the skyrmion crystal is commensurate with the lattice. The different Fermi-surface topology is directly distinguished by an electric probe of angle-resolved photoemission spectroscopy. Furthermore, we show that the SkXs obtained by the variational calculations are also confirmed by numerical simulations on the basis of the kernel polynomial method and the Langevin dynamics for the double exchange model and the simulated annealing for an effective spin model.


I. INTRODUCTION
A wave nature of electrons provides millions of enigmatic problems since the early twentieth century. In particular, an interference among electron wave functions in solids gives rise to unconventional electronic structures and transport properties, which makes considerable advances of modern science technology. For example, specific atomic alignments on lattices bring about unconventional electronic states with massless Dirac fermions, flat-band structures, and quasi-crystal nature as recently observed in the twisted bilayer graphene [1][2][3][4][5]. Another example is found in quantum interference effects through electron dynamics, such as the Anderson localization in disordered systems [6,7], anomalous quantum Hall effect in magnetic topological insulators [8,9], and Josephson effect in superconductors [10,11].
From an energetic point of view, the emergence of these multiple-Q states depends on the model Hamiltonian and the lattice structure [65]. For example, the types of constituent waves in the SkXs are determined by microscopic interactions; the Heisenberg model with the ferromagnetic interaction and the Dzyaloshinskii-Moriya (DM) interaction in chiral (polar) magnets tends to favor the Bloch (Néel) SkX consisting of proper-screw (cycloidal) spirals, while the double exchange (classical Kondo lattice) model on a centrosymmetric square lattice hosts either the Bloch, Néel, or anti-SkX depending on the signs of the anisotropic Ruderman-Kittel-Kasuya-Yosida (RKKY) interactions [66][67][68] in momentum space [69]. Similarly, it is recognized that the phase degree of freedom in the multiple-Q states is closely related to interactions, which brings about the phase transitions between the different types of multiple-Q states [59,70,71]. The long-range chirality interaction leads to a topological transition between the SkX with |n sk | = 2 and the tetra-axial vortex crystal with |n sk | = 0 by changing the relative phase among the constituent waves [59], and magnetic anisotropy arising from the mirror symmetry breaking drives a transition between the SkXs with |n sk | = 1 and |n sk | = 2 by changing the relative phase between xy and z components of the individual spiral [29]. These studies suggest that there is a further possibility to induce a new type of multiple-Q states by considering specific model parameters and lattice structures, which might be useful to stimulate further experimental findings of exotic topological spin crystals beyond the SkX and the hedgehog lattice [65,72]. It is also desirable to provide a way to identify different multiple-Q states experimentally, since it is sometimes difficult to distinguish them by diffraction techniques like the neutron scattering and the resonant x-ray scattering.
In the present study, we investigate the formation of the SkXs in itinerant magnets on a discrete lattice by focusing on a position of the skyrmion core, which is related to the phase degree of freedom in the SkXs. We consider the short-period SkXs in the double exchange (classical Kondo lattice) model with an easy-axis singleion anisotropy on a triangular lattice, where the SkXs with |n sk | = 1 and |n sk | = 2 consisting of the commensurate ordering wave vectors with the crystal lattice are stabilized with and without the magnetic field [26]. By constructing a phase diagram within the variational calculations, we find that the position of the skyrmion core is located at the interstitial site for a small single-ion anisotropy, while that is locked at the onsite for a large single-ion anisotropy. Reflecting the different core positions, the resultant electronic states for the SkXs with |n sk | = 1 are sensitively affected in a different way: the threefold-symmetric Fermi surface appears for a small anisotropy, while the sixfold-symmetric one appears for a large anisotropy. Furthermore, we show that the SkX with |n sk | = 1 is described by superposing of three spiral and sinusoidal waves for a small anisotropy, which is in contrast to the superposition of three spirals for the SkX in chiral/polar magnets. We confirm that a part of the SkXs obtained by the variational calculations are reproduced by the unbiased kernel polynomial method and the Langevin dynamics (KPM-LD) simulations for the large system size [73]. We also perform the simulated annealing for an effective spin model of the double exchange model to discuss the SkX spin texture for a large single-ion anisotropy.
The rest of the paper is organized as follows. After presenting the double exchange model in Sec. II, we show the results by the variational calculations for the double exchange model on the triangular lattice in Sec. III. We discuss the optimal variational parameters and the Fermi surfaces while changing the magnetic field and easy-axis single-ion anisotropy. We examine the stability of the SkXs by using the unbiased KPM-LD simulations in Sec. IV. We also show the result by the simulated annealing for the effective spin model. We summarize our results in Sec. V.

II. DOUBLE EXCHANGE MODEL
We consider an itinerant electron model consisting of noninteracting electrons coupled with localized spins, dubbed the double exchange (classical Kondo lattice) model, on the triangular lattice under the periodic boundary condition. The Hamiltonian is given by where c † iσ (c iσ ) is a creation (annihilation) operator of an itinerant electron at site i with spin σ. The first and second terms represent the kinetic energies of itinerant electrons, where the sums of i, j and i, j are taken over the nearest-and third-neighbor sites on the triangular lattice, respectively. The third term represents the exchange coupling between itinerant electron spins s i = (1/2) σσ c † iσ σ σσ c iσ and localized spins S i , where σ = (σ x , σ y , σ z ) is the vector of Pauli matrices and S i is regarded as a classical spin whose amplitude is normalized as |S i | = 1. The sign of J is irrelevant for classical localized spins. The fourth term represents the chemical potential µ. The fifth term represents the easy-axis single-ion anisotropy (A > 0). The sixth term represents the Zeeman coupling to an external magnetic field along the z direction. We consider the effect of the magnetic anisotropy and the magnetic field only on localized spins for simplicity. Hereafter, we take t 1 = 1 as the energy unit of the double exchange model. We also set the lattice constant a = 1 as the length unit of the triangular lattice.
The other model parameters are chosen so that the SkXs are stabilized in the ground state and their period is commensurate with that of the lattice: t 3 = −0.85, J = 1, and µ = −3.5 [26]. For these parameters, the bare susceptibility of itinerant electrons shows a triple peak structure in momentum space at Q 1 = (Q, 0), with the relatively large ordering vector Q = π/3. The three ordering vectors, Q 1 , Q 2 , and Q 3 , are connected by the threefold rotational symmetry. Consequently, the short period triple-Q states are stabilized in the ground state: the |n sk | = 2 SkX for 0 ≤ H 0.00325, the |n sk | = 1 SkX for 0.00325 H 0.0065, and another triple-Q state with |n sk | = 0 for 0.0065 H in the absence of A [26]. By introducing the easy-axis anisotropy, the stable region for the |n sk | = 1 SkX is extended, e.g., 0.004 H 0.0115 for A = 0.002 [27]. Meanwhile, the detailed multiple-Q nature of the SkXs, especially for the |n sk | = 1 SkX, was not fully shown, as described below. In the following, we discuss the details of the SkXs by explicitly taking into account the phase degree of freedom among the constituent waves with the use of the variational calculations.

III. VARIATIONAL CALCULATIONS
We show the result by the variational calculation in this section. In order to discuss the optimal spin configurations of the SkXs and their modulations while varying the model parameters A and H in Eq. (2), we consider a general form to represent the triple-Q structures, which is constructed by the superposition of the three spin density waves with Q 1 , Q 2 , and Q 3 as where T denotes the transpose of the vector. In the spin configuration, Q νi = Q ν · r i + θ ν (ν = 1-3 and r i is the position vector at site i) and Q νi = Q νi + ψ; θ ν , ψ, a z , andM z are the variational parameters. We choose the specific vorticity and helicity, which are defined by the winding number projected onto the S x -S y plane around the skyrmion core and the relative angle in the S x -S y plane, respectively [15], because they are arbitrary owing to rotational symmetry around the z axis in spin space in the double exchange model in Eq. (2). We do not consider the anisotropic triple-Q state with different amplitudes at Q 1 -Q 3 , which is stabilized for 0.006 H at A = 0 [26], in order to focus on the optimal spin configurations of the SkXs described by Eq. (3) for simplicity.
There are two phase degrees of freedom in Eq. (3). The one is θ ν , which represents the relative phase degree of freedom among the constituent waves. We here consider four sets of (θ 1 , θ 2 , θ 3 ) in the variational calculations: (0,0,0) denoted as the pattern I, (π/3, −π/3, 0) denoted as the pattern II, (π/2, π/2, π/2) denoted as the pattern III, and (5π/6, π/6, π/2) denoted as the pattern IV. The patterns I and II with ν θ ν = 0 represent the SkX spin textures and the patterns III and IV with ν θ ν = 3π/2 represent the meron-antimeron crystal spin textures, as shown in Figs. 1(a)-1(d). The main difference between them is found in the net scalar chirality at zero field (M z = 0). The SkX spin texture exhibits a nonzero net scalar chirality χ sc = R χ R , where R is the position vector for the center of the triangle plaquette, while the meron-antimeron crystal one does not have χ sc owing to the cancellation of positive and negative χ R . It is noted that the patterns III and IV carry a nonzero net scalar chirality forM z = 0 in the external magnetic field, which can also have a quantized skyrmion number [47].
Meanwhile, the difference between the patterns I and II (or patterns III and IV) is found in the position of the skyrmion core defined as the position with S z i = −1. The skyrmion cores in the patterns I and III are lied at the sites, while those in the patterns II and IV are lied at the interstitial sites, as shown in Figs. 1(a)-1(d). In other words, the constitute waves with the same phase (θ 1 = θ 2 = θ 3 ) fix the core at the sites, while those with different phases remove it from the sites.
The other phase degree of freedom in the SkXs is the relative phase between the xy-and z-spin components, ψ, which has been explicitly pointed out in Ref. 29. This phase degree of freedom is related to the nature of the spin density waves for individual wave vectors. The state with ψ = 0 represents the spiral wave [ Fig. 1(g)], while that with ψ = π/2 represents the sinusoidal wave [ Fig. 1(g)] along each Q ν direction. Reflecting the dif- For the above variational states with the patterns I-IV spin textures, we compare the grand potential at zero temperature, Ω = H /N (N is the system size), and determine the lowest energy state. The variational parameter ψ is taken at ψ = πl/240 (l = 0, 1, 2, · · · , 120) for the patterns I and II, while we take ψ = πl/240 (l = 0, 1, 2, · · · , 240) for the patterns III and IV. The other variational parameters, a z andM z , are taken at a z = 0.02l (l = 0, 1, 2, · · · , 300) andM z = 0.002l (l = 0, 1, 2, · · · , 2500). We consider these ordered states for the system size with N = 96 2 and compute Ω by dividing the system into a 144-site unit cell and 64 supercells under the periodic boundary conditions. Figure 2 shows the variational phase diagram of the double exchange model in Eq. (2) while varying A and H. The spin patterns I, II, and IV appear in the phase diagram, while the pattern III is not stabilized in the present parameter region. The dashed curve in Fig. 2 stands for the boundary between the |n sk | = 1 and |n sk | = 2 SkXs. The corresponding variational parameters a z ,M z , and ψ are shown in Figs. 3(a), 3(b), and 3(c), respectively. The behaviors of a z andM z in Figs. 3(a) and 3(b) are intuitively understood from the energy gain in the presence of A and H; a z tends to be larger for large A, since the easy-axis anisotropy aligns the spins along the ±z direction. Similarly,M z tends to be larger for large H, since the magnetic field aligns the spin along the +z direction. Furthermore, we find that ψ also changes depending on H and A, as shown in Fig. 3(c). Here, the type of the spin density wave changes from the sinusoidal wave (ψ = π/2) for a low magnetic field to the spiral wave (ψ = 0) for a high magnetic field. This behavior is attributed to the ψ dependence of the magnetization; the magnetization for ψ = 0 tends to be larger than that for ψ = π/2 in the triple-Q spin texture in Eq. (3) [29].
For the isotropic double exchange model at A = 0, the |n sk | = 2 SkX with the spin pattern II is stabilized for 0 ≤ H 0.00325, while the |n sk | = 1 SkX with the spin pattern IV or II is stabilized for 0.00325 H. The critical field of the transition between SkXs with |n sk | = 1 and |n sk | = 2 is almost consistent with that obtained by the unbiased KPM-LD simulations [26]. It is noted that the |n sk | = 1 SkX for 0.0065 H is metastable, since the anisotropic triple-Q state which is not taken into account in the variational states in Eq. (3) is stabilized. The variational results show that the skyrmion cores in the isotropic double exchange model are locked at the interstitial sites, as shown in Figs. 1(b), 1(d), and 1(f), which indicates sixfold rotational and inversion symmetries of the triangular lattice is broken in both the |n sk | = 1 and |n sk | = 2 SkXs.
Although the |n sk | = 2 SkX in the pattern II region in the phase diagram has a similar spin texture to that in Fig. 1(f), the |n sk | = 1 SkX in the pattern IV region shows a different spin texture from that with ψ = 0 in Fig. 1(d). In this region, the |n sk | = 1 SkX is characterized by ψ ∼ π/3. In other words, the spin texture is described by the superposition of the sinusoidal wave and the spiral wave, which is different from the SkX represented by the superposition of three spirals in chiral/polar magnets [see Figs. 1(a) and 1(b)]. The spin configuration of the pattern IV at A = 0 and H = 0.005 is shown in the left panel of Fig. 3(d), which consists of three types of vortices: one with vorticity −2 around S z +1, another with vorticity +1 around S z +1 and the other with vorticity +1 around S z −1. The corresponding spin chirality distribution is shown in the left panel of Fig. 3(e), where the integration of the topological charge proportional to the spin chirality gives the quantized number of minus one (n sk = −1). Reflecting the breaking of sixfold rotational symmetry, the Fermi surface in the ordered state is threefold rotational symmetric, as shown in the left panel of Fig. 3(f).
By introducing A, the optimal spin configuration of the |n sk | = 1 SkX in the pattern IV region is mainly replaced by the pattern II spin configuration while keeping |n sk | = 1 except for the narrow pattern I region sandwiched between the patterns IV and II regions. The spin texture in the pattern II is characterized by small ψ, as shown in Fig. 3(c), whose spin and chirality configurations are shown in real space in the middle left panels of Figs. 3(d) and 3(e), respectively. Compared to the left and middle left panels in Figs. 3(d) and 3(e), one can find that both the spin and scalar chirality textures seem to be similar with each other. Thus, it is difficult to distinguish the patterns II and IV in the external magnetic field. Meanwhile, the important observation is that the skyrmion core is still placed at the interstitial site when introducing A, although the pattern I appears in the narrow region between the patterns II and IV. As the skyrmion core is located at the interstitial site, the threefold-symmetric Fermi surface is also obtained in the pattern II, as shown in the middle left panel of Fig. 3(f).
While further increasing A, the pattern II state turns into the pattern I state for both the |n sk | = 1 and |n sk | = 2 SkXs, as shown in Fig. 2. In the region for the |n sk | = 1 SkX, the spin texture is characterized by ψ = 0, which indicates that the individual spin density wave is described by the spiral state. In this case, the skyrmion core is lied at the onsite position, as shown in the middle right panel of Fig. 3(d). Accordingly, the scalar chirality is also distributed in a sixfold-symmetric way, as shown in the middle right panel of Fig. 3(e). In the end, the Fermi surface is sixfold symmetric in the middle right panel of Fig. 3(f) [74]. This sixfold Fermi surface is again broken by increasing the magnetic field H, where the pattern I state changes into the pattern II one. The skyrmion core in the pattern II state is located at the interstitial sites as shown in the right panels of Figs. 3(d) and 3(e), and then the sixfold rotational symmetry of the Fermi surface is slightly broken in the right panel of Fig. 3(f).
For low-field H, the optimal spin configuration of the |n sk | = 2 SkX in the pattern II region is almost the same against A up to A 0.0026, as shown in Fig. 2. Similar to the |n sk | = 1 SkX, the skyrmion core position moves to the onsite position for large A, which is denoted as the pattern I. However, the spin texture in the |n sk | = 2 SkX in Eq. (3) has higher energy than the other spin ansatz in the form of S i ∝ (cos Q 1i , cos Q 2i , a z cos Q 3i ) with |n sk | = 2 when the unbiased KPM-LD simulations are performed [27]. This is distinct from the situation where the |n sk | = 1 SkX in Eq. (3) is stabilized under the magnetic field, as discussed above; indeed, we show that the |n sk | = 1 SkX expected from the variational calculations remains stable in the presence of A even after carrying out the KPM-LD simulations, as will be discussed in the next section [27].
The above results indicate that the skyrmion cores are located at the interstitial sites without the anisotropy, while they move to the onsite when increasing A. This behavior is reasonable because the easy-axis anisotropy favors the large z-spin component, whose energy gain is maximized when the skyrmion core with S z i −1 is placed at the site. Although it is usually difficult to control A by external stimuli, one can also change the position of the skyrmion core by applying the external magnetic field for materials with the large easy-axis anisotropy. The change of the core position can be directly detected in k-resolved angle-resolved photoemission spectroscopy as a consequence of different symmetry of the Fermi surface. It is noteworthy such information to identify complicated spin textures by an electric probe has been achieved by the noncollinear magnetoresistance [75,76] and the spectroscopic imaging scanning tunneling microscopy measurement [77], since the formation of the (multiple) spin density waves affects local electronic states and causes charge density waves [78]. Our results indicate that angle-resolved photoemission spectroscopy becomes another electric probe of the SkX spin textures in a complement way. Moreover, as the antisymmetric band structure with respect to k is related with the nonreciprocal transport even without the spin-orbit coupling [79][80][81][82], our results indicate that the control of the transport properties is possible in the SkXs in itinerant magnets.

IV. LANGEVIN DYNAMICS SIMULATION
To investigate the stability tendency of the |n sk | = 1 SkX in the previous section, we perform the KPM-LD simulation, which is an unbiased numerical simulation based on Langevin dynamics [73] with the kernel polynomial method [83]. As this method has an advantage of performing a large-scale calculation compared to the standard Monte Carlo simulation combined with the direct diagonalization, it has been applied to similar models consisting of free fermions and one-body boson field [26,61,73,[84][85][86][87][88][89][90]. The simulation is performed at zero temperature for the system size with N = 96 2 under the periodic boundary conditions. In the kernel polynomial method, we expand the density of states by up to 2000th order of Chebyshev polynomials with 16 2 random vectors [91]. In the Langevin dynamics, we use a projected Heun scheme [92] for 1000-3000 steps with the time interval ∆τ = 2. Figure 4(a) shows the real-space spin configuration at A = 0 and H = 0.006 obtained by the KPM-LD simulations. The spin texture well corresponds to that obtained by the variational calculations in the left panel of Fig. 3(d): it consists of three types of vortices, one with vorticity −2 around S z +1, another with vorticity +1 around S z +1 and the other with vorticity +1 around S z −1. We also reproduce the spin configuration of the pattern II in the middle left panel of Fig. 3(d) by performing the KPM-LD simulations at A = 0.002 and H = 0.005 in Fig. 4(b). The scalar chiralities in real space in Figs. 4(d) and 4(e) correspond to those in left and middle left panels of Fig. 3(e) as well. Similarly, the n sk = 2 SkX at A = 0 is obtained by the KPM-LD simulations [26].
However, we could not obtain the spin texture corresponding to that in the middle right panel of Fig. 3(d) in the KPM-LD simulations, since the simulations are easily trapped into the metastable state with the multidomain configurations in terms of the spin and chirality, which makes the core position ambiguity. This is presumably attributed to the small energy scale under the large single-ion anisotropy. Instead of that, we analyze an effective spin model derived from the double exchange model, where the first four terms in Eq. (2) are replaced by the bilinear RKKY interactionJ and the positive biquadratic interaction K (see Ref. 93 for the model in details). By taking the model parameters as K/J = 0.5, A/J = 0.2, and H/J = 0.4, adopting the same ordering vectors, and performing the simulated annealing following the manner in Ref. 93, we confirmed that the skyrmion core is located at the site for a relatively large A, as found in the middle right panel of Fig. 3(d). The real-space spin and scalar chirality configurations obtained by the simulated annealing are shown in Figs. 4(c) and 4(f), respectively.

V. SUMMARY
To summarize, we have investigated the optimal spin configurations of the SkXs in itinerant magnets by tak-ing into account the phase degrees of freedom among the constituent waves. We find that the core positions of the SkX consisting of the commensurate ordering wave vectors with the crystal lattice are placed at the interstitial sites for small magnetic anisotropy, while those are found at the sites for large magnetic anisotropy by the variational calculations for the double exchange model on the discrete triangular lattice. We also find that the SkX with |n sk | = 1 consists of the superposition of the triple sinusoidal and spiral waves for small magnetic anisotropy, while that is described by the triple spiral waves for large magnetic anisotropy. We discuss the possibility of the control of the nonreciprocal transport by applying the magnetic field (or changing the single-ion anisotropy), since the core positions at the interstitial site (onsite) break (hold) the sixfold rotational symmetry, leading to the threefold(sixfold)-symmetric band deformation. On the basis of the different Fermi-surface geometry resulting from the different core position, we propose that angle-resolved photoemission spectroscopy can be a useful electric probe to obtain information about the real-space spin textures. Finally, we confirm that the spin textures expected from the variational calculations are obtained by the unbiased KPM-LD simulations for the double exchange model or simulated annealing for the effective spin model.
The present results provide information to examine the detailed symmetry of not only the SkXs and but also the hedgehog lattices with the commensurate magnetic modulations where the discreteness of the lattice structure is important [59,70]. Especially, the SkXs and the hedgehog lattices stabilized by the spin-charge entanglement in itinerant magnets might be promising, since the commensurate locking often happens owing to the local gap formation of the electronic band structure, as found in the chiral soliton lattice [94,95]. In addition, the results would be useful to investigate the effect of the intrinsic lattice pinning in the current-driven dynamics [96,97].